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... A simple two-layer variable infiltration .capacity (VIC-2L) land surface 
for incorporation in general circulation models (GCMs) is 
described. The model Consists of a two-layer characterization of the soil within a 
ULwl grid cell, and uses an aerodynamic representation of latent and sensible 
heat fluxes at the land surface. The effects of GCM spatial subgrid variability of 
soil moisture and a hydrologicaily realistic runoff mechanism are represented in 
the soil layers. In the upper layer, the spatial distribution of infiltration and soil 
moisture capacities is included. The lower layer is lumped spatially and uses a 
nonlinear drainage representation. The model partitions the <£ea of interest into 
multiple land surface cover types; for each land cover type the fraction of plant 
roots in the upper and lower zone is specified. Evaporation occurs via caLpy 
evaporation, evaporation from bare soil, and transpiration, which is represented 
using a canopy and architectural resistance formulation. The model was tested 
using long-term hydrologic and climatalogical data for Kings Creek, Kansas to 
estimate and validate the hydrological parameters. Surface flux data from three 
Satellite Land Surface Climatology Project Field Experiment 
(Mb-,) intensive field campaigns in the summer and fall of 1987 in central 

^ 0I § An Sfo- B *^an Amazonian Climate Observation Study 
( ABRACOS) in Brazil were used to validate the model-simulated surface energy 
fluxes and surface temperature. 


In addition, a derived distribution approach which accounts for the effects of 
subgrid scale spatial variabilities of precipitation on surface energy fluxes, soil 
moisture, and runoff production was developed for an extended version of VIC- 
2 T i , n } s 0 jf - *“ e derived distribution approach differs from pixel-based approaches 
which discretize precipitation over a spatial domain, and from previous statistical 
approaches that combine the point precipitation distribution with the point 
statistical distribution of Selected land Surface characteristics. The results of the 
denved distribution method are compared with those obtained using an 
exhaustive pixel-based approach, and the results obtained by applying uniform 
spatiafly averaged precipitation to the VJC-2L model. Under most conditions, 
the denvew distnbution approach gives good approximations to the pixel-based, 
approach, and is superior to the constant precipitation approach for surface 
fluxes, surface temperature, runoff, and soil moisture. Finally, VIC-2L sensitivity 
of predictions to model parameters Were explored for two different climate 
regimes using both fractional factorial and one-at-a-time sensitivity analyses. 
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CHAPTER 1 INTRODUCTION 


1.1. Motivation for and dissertation structure 

The parameterization of the hydrologic and thermal characteristics of the 
land surface is important for general circulation models (GCMs) used fc r cl mate 
prediction and weather forecasting (e.g., Dickinson 1991, Dickinson and 
Kennedy 1991, Wood 1991, Shuttleworth 1991a, Henderson-Sellers and Pitman 
1992, among others). Garratt (1993) reviewed the land surface and boundary- 
layer treatments in some 20 GCMs through sensitivity studies of climate 
simulations, and found that the regional and global climate are most sensitive to 
albedo, surface roughness length, vegetation coverage, and Soil moisture 
distribution. The. inclusion of a canopy scheme which allows more reasonable 
consideration of effects of albedo, roughness, and soil moisture significantly 
improved the simulated climate. It . also facilitated studies of the effect of 
deforestation on climate (see, for example, studies of the regional impact of 
Amazonian deforestation by Sud et al. (1990), Pitman et al. .(1993), Henderson- 
Sellers et al. (1993), and Eltahir and Bras (1993), among others). 

Although complicated canopy schemes such as BATS (Biosphere- 
Atmosphere Transfer Scheme, Dickinson et al. 1986), SiB (Simple Biosphere 
model, Sellers et al. 1986), and BEST (Bare Essentials of Surface Transfer, 
Pitman et al. 1991), and simpler . canopy models (e.g., Noilhan and Planton 
1989) have been implementated into GCMs to give more reasonable climate 
simulations, precipitation, evaporation, soil moisture, and surface temperature 
axe still not, in general, well simulated. Examination of systematic errors in 
GCM climate simulations, particularly over land areas, is an active research 
area. Verseghy et al. (1993) found that the assumption made in many land 
surface schemes that excess surface water is removed immediately from the land 
surface system can result in .substantial overestimates of surface temperatures in 
continental interiors. Garratt (1993), in his review of sensitivity studies, 
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suggested that there is a need for a more realistic land surface representation in 
GCMs, particularly with respect to (1) surface fluxes at the appropriate 
horizontal resolution, (2) surface runoff and canopy interception processes, and 
(3) the spatial distribution of rainfall within a GCM grid. 

In this dissertation, a two-layer variable infiltration capacity (VIC-2L) 
model that attempts to incorporate the above three features suggested by 
Garratt (1993) is developed. The model includes a canopy layer, the effects of 
spatial Subgrid variability of soil moisture with a hydrologically reasonable runoff 
mechanism, and the influence of the subgrid spatial distribution of rainfall. The 
development of the VIC-2L model is described in Chapter 2, and is evaluated 
using observed data from the FIFE (central Kansas, USA) and ABRACOS 
(Amazonia) sites in Chapter 3. These two sites are sufficiently small that the 
assumption of constant precipitation within the measurement area is reasonable* 
Chapter 4 explores the sensitivities of the VIC-2L model parameters. In Chapter 
5, an extension of the VIC-2L model with a one-diinensional derived distribution 
representation of spatial subgrid variability in precipitation is described. The 
results for versions of the VIC-2L model with constant and Spatially varying 
precipitation are compared and evaluated. Conclusions and recommendations of 
future work are given in Chapter 6. In the remainder of this chapter, several 
land surface parameterizations developed previously are reviewed. Land surface 
schemes that include the spatial variability of rainfall within a grid area are 
reviewed in Chapter 5. 

1.2. Background 

The problem of how to represent land Surface processes in general 
circulation models that are used for climate simulation and numerical weather 
prediction has drawn the interest of climate modelers, and increasingly, 
hydrologists and systems ecologists. Early generation GCMs did not include 
representations of land surface hydrology, instead they used prescribed surface 
wetness and temperature, and thus could not account for the feedbacks between 
the land surface and atmosphere. In retrospect, such representations have 
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proved useful only for examining the limiting cases: perpetually wet and dry 
surfaces (Shukla and Mintz, 1982). These early studies confirmed the 
importance of including the interactions between the land surface hydrology and 
atmosphere. 

Notwithstanding the desirability of better representing the land surface in 
GCMs, the spatial scale of the GCM ’’grid box” is so large (typically at least 100 
km) that only relatively simple models can be justified, especially considering 
the other sources of uncertainty in climate and weather prediction (e.g., cloud 
physics, ocean circulation). Manabe et al. (1969) followed this logic in using • 
Budyko’s ’bucket’ model to represent the land surface hydrology at the global 
scale. The bucket model assumes that all rainfall, is infiltrated until soil moisture 
capacity is exceeded, whereafter . the excess precipitation becomes runoff. 
Manabe ’s. bucket model uses a simple relationship between actual and potential 
evaporation which is often referred to as a beta function. Milly (1992) pointed 
out the conceptual inconsistency in calculating the potential evaporation and the 
coefficient ’’beta” for this kind of model formulation by using the modeled 
surface temperature to evaluate potential evaporation. He argued that the 
appropriate temperature to use for evaluating potential evaporation is that of a 
freely evaporating surface, and described, two approaches to remove the. 
inconsistency. The bucket model is clearly simplistic with respect to infil tration 
and runoff production, in addition to evaporation, and ignores vegetation effects 
on evapotranSpiration, except to the extent that the beta function, acts as a 
surrogate. Another problem with the bucket model, is that in its most common 
implementation, the parameters are assumed to be constant over the globe, 
although this assumption can be relaxed using suitable d ata, such as globally 
varying soil water holding capacities. 

While the shortcomings of the bucket model and related simple 
representations of the land surface, such as the two-layer surface model used in 
the GISS GCM (Hansen et id. 1983) are clear, the pathway toward more 
realistic parameterizations is less obvious. The difficulty of estimating 
parameters globally, as well as the desire to keep the complexity of the land 
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surface representation compatible with that used to represent other elements of 
the atmosphere-land-ocean system, are important considerations in determining 
what form the next generation of land surface hydrology models for GCMs 
should take. Mintz (1984) and Rowntree (1988) reviewed the impact of land 
surface boundary conditions on . simulated climate and pointed out that the 
atmosphere is sensitive to land surface evapotranspiration, which is largely 
affected by changes in available soil moisture or in albedo.. Their review 
stimulated the incorporation of improved of canopy and soil formulations in 
land-surface schemes in GCMs (Garratt 1993). 

One direction that has been pursued is to improve the representation of 
soil moisture dynamics, and especially to represent vegetation interactions with 
the soil column and the atmosphere explicitly so that effects of biosphere-climate 
interactions can be studied. The resulting class of models is known as soil 
vegetation atmosphere transfer schemes (SVATS). Among the SVATS that 
have been developed for GCM use are BATS (Dickinson et al. 1986, Dickinson 
et al. 1993) and SiB (Sellers et al. 1986). A distinguishing feature of SVATS, 
which is evident in both BATS and SiB, is that they have a high level of 
vertical resolution and structure, but a low level of horizontal resolution (Wood 
1991). For instance, the parameters for the soil and vegetation properties are 
assumed constant within, a GCM gnd, thus ignoring spatial heterogeneity. In 
addition, most SVATS use a "flat earth” representation of the land Surface 
which neglects, the effects of topography on runoff production and soil moisture 
dynamics. Because of their, emphasis on vertical structure, SVATS arguably do 
a better job of partitioning incoming solar radiation into latent and sensible . heat 
than they do in accounting for soil moisture dynamics and runoff production.. 

An alternative line of investigation is to develop simpler land surface 
models that still incorporate important features of the governing hydrological 
processes. For example, Xue et al. (1991) simplified SiB in three aspects. 
These are the parameterization of the diurnal variation of surface albedo, the 
effect of root zone soil moisture on stomatal resistance, and the surface stress 
and fluxes of heat and moisture between the top of the vegetation can opy and 
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the atmospheric reference level. -After these simplifications, the two vegetation 
stories in SiB become oiie layer. With negligiblelosS of accuracy, the number of 
parameters of SiB model Was reduced by Xue et al. (1991) from 44 to 21. 

Mahrt and Pan (1984) developed a two-layer soil hydrological scheme for 
use in GCMs and numerical weather prediction models. There are three major 
features in their scheme. First, it has a thin upper layer used to. represent the 
large vertical moisture gradient near the surface during evaporation. Second, it 
includes a strong dependence of hydraulic diffusivity on the vertical distribution 
of soil moisture. Finally, it estimates surface evaporation by using near-surface 
soil water flux information. They also suggests that transpiration can be 
estimated by relating it to layer^averaged moisture and- potential evaporation if 
Vegetation is present. 

Abramopoulos et al. (1988) developed a. simple land surface scheme that 
includes multiple soil layers with specified hydraulic conductivity and matric 
potential functions rather than a constant diffusivity to describe the soil water 
dynamics by using Darcy’s law. Vegetative resistance, evaporation from 
intercepted precipitation and dew, evaporation from bare soil, and transpiration 
are all explicitly represented. Abramopoulos et al. (1988) compared the 

evaporation predicted using am areally Weighted average of a heterogeneous land 
surface with the evaporation obtained using area- Weighted vegetation and soil 
characteristic parameters, and found that it is better to average the subgrid 
fluxes than to average the soiLand vegetation parameters. 

Noilhan . and Pianton (1989) developed a simple land surface 
parameterization for meteorological models. Unlike BATS and SiB where the 
canopy and soil temperatures can vary, they treated the system as isothermal 
with equal canopy and the upper soil layer temperatures. Also, their model 
does not calculate soil fluxes beneath the vegetation. BATS and SiB both 
calculate detailed fluxes from the soil beneath the canopy, the open areas 
between the canopy, and the canopy itself. 

Pan (1990) used the bucket concept to estimate the actual evaporation 
with vegetation lumped with soil in the National Meteorological Center medium- 
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range forecast (NMC MRF) model. He overcame the inconsistency problem 
noted by Milly (1992) by calculating the potential evaporation with a 
hypothetical temperature from a wet surface using the Mahrt and Ek (1984) 
method, rather than with the unsaturated surface temperature when the soil, is 
dry. When vegetation is present, potential evaporation is calculated using 
Monteith’s minimum resistance concept. 

Pitman et al. (1991) and Yang and Pitman (199u, developed a land 
surface scheme which uses simplified (as compared with BATS) albedo and 
stomatal resistance formulations, and explicit representation of frozen soils. In 
addition, it calculates infiltration, runoff, and soil evaporation following 
Eagleson’s (1970) approach. The canopy iayer is formulated as a nonisothermal 
system like BATS, although the algorithm for. canopy temperature is different. 

Sielert et al. (1992) developed a soil-vegetation model for use in a 
mesoscale atmospheric model. The canopy model is based on the work of 
Deardorff (19.78) and Dickinson (1984). The soil model is based on Sievers et al. . 
(1983) which includes a complicated treatment of heat and moisture transport 
within the Soil. The soil temperature and moisture content are represented by 
two coupled differential equations^ 

Pollard and Thompson (Bonan et al. 1992, Bonan et al. 1993) developed 
a six-layer soil model for GCMs that explicitly represents the effects of Soil 
freezing and thawing on latent heat flux. The snow cover is represented . by 
three-layers. The canopy layer is formulated in a similar way to BATS and SiB, 
but is less complex. The scheme has. been used to study the effects cf boreal 
forests on climate in the National Center for Atmospheric Research co mmuni ty 
climate model CCM1 (Bonan et al. 1992). 

Hondo and Watanabe (1992) developed a multilayer energy budget model 
for a vegetation canopy. The model represents the energy budget of leaf 
surfaces, the ground surface, and turbulent and radiative transfer processes 
within the canopy. The vegetation is partitioned into 50 layers and the energy 
budget equation for each layer is solved for steady and horizontally homogeneous 
flow. 
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Numerical weather, modelers have begun to include simplified 
parameterisations of iand surface processes into operational numerical Weather 
prediction models. Examples are . the models of the European Center for 
Medium-range Weather Forecasting (ECMWF), the U.S> National 
Meteorological Service (NMC), the UK Meteorological Office, the French 
Meteorological Service (Direction de la Meteorologie), and the Japan 
Meteorological Agency (Blondin 1991). These models have many features in 
common. They all simulate the diurnal cycle, have their first atmospheric level 
inside the constant flux layer (CFL), use Monin-Obukhov similarity theory, 
and use only one roughness length for momentum, heat, and moisture. The 
surface moisture flux, however, is treated differently in each model. 

Although the models discussed above have various degrees of complexity 
in their canopy and soil representations, none of them considers subgrid scale 
spatial variability in either meteorological inputs to the. land surface or land 
surface characteristics. Representation of heterogeneities in terrain, soil, 
Vegetation, and precipitation at scales smaller than those resolved by GCMs has 
been a relatively recent concern. Warrilow et al. (1986) took account of spatial 
variations in rainfall by assuming that only a fraction of a grid cell receives 
rainfall, and that Within the fraction covered, rainfall depth is distributed 
exponentially. However, in their parameterization, the interception of rainfall 
was assumed to be uniformly distributed over the entire grid cell, and thus 
evaporation was overestimated (Lean and WaxriloW 1989). Following .the 
formulation of Dolman and Gregory . (1992), Lean and Rowntree (1993) 
incorporated a new interception representation Which assumed that interception 
occurred only over a fraction of . the grid cell within which the rainfall Was 
governed by an exponential distribution. 

Wetzel and Chang (1988) incorporated the effect of subgrid variability in 
soil moisture into the evaporation process ~by using a statistical distribution to 
represent the soil moisture within a grid cell. The .soil column in their model 
consists of three layers With a thin surface layer and two thicker sublayers. The 
first sublayer is assumed to contain 50% of all plant roots and is used to 
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represent the diurnal Variation of soil moisture; and the second sublayer 
Contains the other half of the roots and is used to represent time scales of one 
day or longer. Evaporation is evaluated based on the Ohm’s law analog form. 
The grid-cell-average evapotranspiration was then summed based On weighting 
factors obtained. from the statistical distribution Of the soil moisture. 

Avissar and Pielke (1989) investigated the subgrid scale variability 
associated with land surface heterogeneities in a mesoscale model by classify ing 
the surface into similar homogeneous patches. After regrouping into Subgrid 
classes, a parameterization for homogeneous surfaces was applied. The total- 
fluxes of energy of each grid cell were then evaluated according to the 
distribution of the different subgrid types within the grid cell. Their Work 
showed that spatial heterogeneity in vegetation can have significant effects on 
temperature and precipitation* 

Entekhabi and Eagleson (1989) prescribed the subgrid spatial Variability 
of soil moisture and storm precipitation statistically and derived expressions for 
hydrologic fluxes based on assumed subgrid soil and precipitation Variability. 
However, their analysis is limited to specific assumed statistical distributions. 
The scheme was implementated into the NASA Goddard Institute for Space 
Studies . (GISS) GCM by Johnson et al. (1993). 

Pielke et al* (1991) illustrated the range of observed spatial variabilities of 
landscape characteristics using observational evidence from field and satellite 
data. They investigated the effects of spatial variability in land surface 
characteristics on lower tropospheric fluxes of energy in the absence of clouds. . 
On this basis, they argued that the influence of mesoscale landscape spatial 
variability on the atmosphere should be explicitly parameterized. 

Avissar (1991) numerically aggregated grid scale surface fluxes for nine 
classes of stomatal resistance. In the aggregation, two stages of scaling were 
used. The first was from the leaf scale to the patch scale; the second was from 
the patch scale to the GCM scale. 

Avissar (1992) used a statistical-dynamical approach to investigate 
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subgrid Scale heterogeneity in stomatal resistance, where a. probability density 
function rather than a single representative value was used. Each term in the 
land surface energy budget was expressed as an integral of a probability density 
function of stomatal resistance. Five different probability density functions for 
stomatal resistance were explored. Comparisons between the results obtained by 
the statistical-dynamical approach and by use of a single representative stomatal 
resistance showed that there were large absolute and relative differences due to 
the nonlinearity of ! and- atmosphere interactions. The observations of stomatal 
conductances in a potato field in New Jersey during the summer of 1989 
supported use of a two parameter lognormal distribution to describe the 
distribution of stomatal conductances. Bonan et al. (1993) also used the 
statistical-dynamic approach to study the effects of subgrid scale heterogeneity in 
leaf area index, minimum and maximum stomatal resistance's, _and soil moisture 
on grid scale fluxes. 

Koster and Suarez (1992a) proposed a mosaic strategy, similar to that 
used by Avissar and Pielke (1989), to account for the effects of different 
vegetation on surface energy fluxes. In the mosaic strategy, each vegetation 
type is represented by a tile, which is coupled independently to the atmosphere. 
The rest of the canopy and soil parameterization is similar to SiB (Sellers et al. 
1986). Koster and Suarez (1992b) extended their earlier mosaic strategy by 
separating the total turbulent flux into latent and sensible heat components. In 
addition, they compared the mosaic approach with a mixture approach in which 
different. Vegetation types were assumed to be homogeneously mixed over a GCM 
grid cell. 

The land surface scheme for the Canadian Climate Center GCMs was 
developed by Verseghy (1991) and Verseghy et al. (1993). In this model, three 
soil layers are used to represent both thermal and moisture regimes. When snow . 
is present, it is assumed that the entire area is covered if the snow depth is. 
greater than an assumed surface roughness height of 0.10 m; otherwise, only a 
fraction of ground area is covered by snow with a fixed depth of 0.10 in.. At the 
end of each time step, the temperature and moisture content are averaged by 
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the relative fraction of snow covered area. In addition, the model calculates the 
canopy and ground surface temperatures separately. As in the strategy adopted 
by Avissar and Pielke (1989), each grid cell is divided into four, different 
subareas. They are., bare soil, snow-covered, VegetatiOn-cOvered, and 
vegetation-and-snow covered areas. Within the vegetation-covered area, average 
canopy parameters based on the presence of four different vegetation types are 
calculated. 

Ducoudre et al. (1993) developed a set of parameterizations of the 
hydrologic exchanges at the land/atmosphere interface within a GCM. The * 
model allows seven different vegetation classes to be present simultaneously 
within the 3ame grid cell. The total latent heat flux transferred to the 
atmosphf *e is obtained by taking the average of evaporation from bare soil, 
transpirc-ion and interception loss from each of the 7 vegetation types. In the 
soil parameterization of the model, rain fills the soil column from top to bottom, 
and water is removed from the Closest level where it is available. 

Famiglietti and Wood (1994a) developed a local water and energy balance 
model which is appropriate for a stream catchment, but could be generalized to 
a region such as a GCM grid cell. The model partitions the land surface into 
bare Soil, wet Canopy, and dry canopy using a pixel-based representation Of the 
land Surface derived from digital elevation data. Spatial variability is explicitly 
incorporated by discrete variation of the model parameters and inputs over the 
Spatial domain. The local fluxes of each, grid element are aggregated either _ 
explicitly or by statistically aggregating the local fluxes through integration over 
their respective spatial probability density functions (Famiglietti and Wood 
1994b). 

One of the major complications in developing and testing land surface 
parameterizations for GCMs is that validation opportunities are few. A recent 
exception is the work of Betts et al. (1993) who compared surface energy fluxes 
and soil temperatures predicted by the ECMWF land Surface model with data 
collected at the FIFE site in central Kansas du ring the summer of 1987. Their 
work identified problems with four components of the model: the incoming 
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shortwave radiation for clear sky conditions, the ground heat flux, surface 
evaporation, and the -entrainment at the top of the boundary layer. They 
concluded that the relatively thin (70 mm) top ground layer in the model was 
the main reason for the errors in ground heat flux, and that deficiencies in the 
surface evaporation algorithm were responsible for large errors in the latent heat 
flux. 

1.3. Research objectives - 

The objectives of this research are: (1) to develop a simple land surface 
scheme appropriate for GCMs that represents spatial variability in soil 
characteristics, vegetation, and precipitation, and simulates explicitly direct 
surface runoff and subsurface runoff; (2) to evaluate the model using observed 
data; and (3) to explore the model parameters using sensitivity analysis. 

The scheme to be developed is a generalization of the VIC model 
described by Wood et al. (1992) and implemented in the GFDL-GCM by Stamm 
et al. (1994). The new model consists of. a simple two-layer characterization of 
the soil Column, and uses an aerodynamic representation of the latent and 
sensible heat fluxes at the land Surface. The soil moisture algorithm is a 
generalization of the Arno model (Francini and Pacciani 1991) in which the 
infiltration, evaporation, soil moisture, and runoff generation vary within an 
area (or within a grid cell in GCMs). The infiltration algorithm in the VIC 
model can be interpreted within the context of a spatial distribution of soils of 
varying infiltration . capacities. It allows, different , types of vegetation to be 
present simultaneously. In addition, it accounts for the spatial variability in 
precipitation. Simplifications of the Arno model using the traditional beta 
function representation of evapotranspirationj have previously been incorporated 
in the Geophysical Fluid Dynamics Laboratory (GFDL) GCMs by Stamm et al. 
(1994), and in the Max Planck Institut GCM by Dumenil and Todini (1992). 

There are major differences between the two-layer VIC model .to be 
described here and the earlier versions incorporated in the GFDL and MPI 
GCMs. The most important differences are the following: (1) both of the earlier 
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schemes have a single soil layer, and neither explicitly represents vegetation in 
the surface energy flux; (2) both of the earlier models distribute precipitation 
uniformly over the grid cell, Stamm. et al. (1994) concluded that ”... the results 
over North American and Eurasia [suggest] the need to represent the surface 
hydrology with a two layer soil system ...”. 


CHAPTER 2 DESCRIPTION OF TWO-LAYER VIC MODEL 


In this chapter, the structure and the characteristics of the two-layer 
variable infiltration capacity model (Liang et- al. 1994) axe described. The 
definition of and dimension of each variable used in this and the following 
chapters axe given in the list of symbols at the start of this dissertation. 

2.1. Introduction 

The model developed here characterizes the subsurface as consisting of 
two soil layers. The Surface is described by N+l land cover types, where n = 1, 
2, * • • , N. represents N different types of vegetation, and n = N+l represents 
bare soil. There , is no restriction on the number of vegetation types in the 
model, but in the interest of model parsimony, N will almost always be less 
than 10. The vertical and. horizontal characterizations are shown schematically 
in Fig. 2.1. The land cover types axe specified by their leaf axea index (LAI), 
canopy resistance, and relative fraction of roots in each of the two soil layers. 
The evapotranspiration from each vegetation type is characterized by potential 
evapotranspiration, together, with, canopy resistance, aerodynamic resistance to 
the transfer of Water, and architectural resistance. Associated With each land 
cover class is a single canopy layer, soil layer 1 (upper zone) and soil layer 2 
(lower zone) . 

The upper layer (soil layer 1) is designed to represent the dynamic 
behavior of the soil that responds to rainfall events, and the lower layer (soil 
layer 2) is used to characterize the seasonal soil moisture behavior. The lower 
layer only responds to rainfall when the upper layer is wet, and thus cm 
separate the subsurface flow from storm quick response. Roots, can extend to 
layer 1 or layers 1 and 2, depending on the vegetation and soil type. For the 
bare soil class, there is no canopy layer. In the present form of ,the model, the 
soiL-characteristics (that is, the distribution of water holding capacities, as 
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described below) are the same for all laid cover classes. However, each cover 
class may have different soil moisture distributions at each time step. 
Infiltration, -drainage of moisture from layer 1 to layer 2, surface, runoff, and 
subsurface runoff are computed for each cover type. The total latent heat flux 
transferred to the atmosphere, total sensible heat and ground heat fluxes, the 
effective surface temperature, and the total surface runoff and subsurface runoff 
are then obtained by summing, over all of the surface cover classes. 

2.2. Evapotranspiration 

Three types of land-atmosphere vapor transport axe considered in the 
model. They are evaporation from the canopy layer of each vegetation class, 
transpiration from each of the vegetation classes, and evaporation from bare soil. 
Total evapotranspiration over a grid, cell (or an area) is. computed as the sum of. 
the canopy, vegetation, and bare soil Components, weighted by the respective 
surface cover area fractions (see Section 2.6). 

^ The maximum canopy evaporation rate for the nth Surface cover class, 
E c [n] , is specified as 



0 


w iM 


W: 


imi 


r wN+r 0 [n" 


( 2 . 1 ) 


In Eq. (2.1), the argument n refers to the vegetation surface cover class index; 
throughout the remainder of this thesis the dependence of many of the surface 
and subsurface characteristics on surface cover class is implied by this argument 
even if not noted specifically. In Eq. (2.1), Wj[n] is the amount of intercepted 
water in storage in the canopy layer, Wj m [n] is the maximum amount of water 
that the canopy layer *.an intercept, Ep[n] is the potential evaporation rate from 
a thin free water surface (Shuttlewofth 1993), r Q [n] is the architectural 
resistance that is due to the variation of the gradient of specific humidity 
between the leaves and the overlying air in the canopy layer (Saugier and Katerji 
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1991), and r w [n] is the aerodynamic resistance to the transfer of Water. The 
power of 2/3 in Eq. (2.1) is used according to Deardorff (1978). The form of Eq. 
(2.1) is a ”beta” representation* where Ep[n] can be obtained by either 
Penman’s formulation (1948) or Penman-Monteith’s formulation for a free water 
surface (Shuttleworth 1993). Penman’s formulation Can be expressed as, 


* r_i _ AR nW + Pa c p( e S- e )/ r a[ n l 

E P [nJ - . W -Le-(A + T) 


( 2 . 2 ) 


where A is the rate of change of saturation vapor pressure with temperature, 
R n [n] is the net radiation, p & is the mass density of air, Cp is the Specific heat of 
air at constant pressure, e s and e are the saturated vapor pressure and vapor 
pressure respectively, r a [n] is the aerodynamic resistance to the momentum 
transfer in the atmosphere (subscripts ”a” could beJ’w’’, ”h’’, etc.), p w is the 
density of liquid Water, L e is the latent heat of vaporization, and 7 is the 
psychr metric constant. The Penman-Monteith’s formulation for a free water 
surface is expressed similarly to Eq. (2.2), except that the available-energy is 
substituted for net radiation (Monteith and Unsworth 1990). 

The maximum amount of water intercepted by the canopy can be 
calculated using the formulation of Dickinson (1984), 

W im M = K L *LAI[n,m] (2.3) 

where is a constant, taken to be 0.2 mm following Dickinson (1984), and 
LAI[n,m] is leaf area index for the nth surface cover class in month m. The leaf 
area index of a canopy is the projected leaf area pet unit ground surface area 
(Campbell 1977). 

The aerodynamic resistance to the transfer of water r w [n] (he., subscript 
”a” becomes ”w”) is calculated as (Monteith and Unsworth 1990): 
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IwW - -ey.nl U D ( Z2 ) (14) 

where u n (z 2 ) is the wind speed over the nth Surface cover class at level z 2 [n], 
and C w [n] is the transfer coefficient for water which is estimated taking into 
account atmospheric stability (Louis 1979) as follows: 

Cw[nl - 1.351 x a2[n] x F w [n] (2.5) - 


where 



. K 2 

z 2 l n i ~ d ofr 
• V 


)] 2 


( 2 . 6 ) 


is .the drag coefficient for the case of. near-neutral stability, K is von Karman’s . 
constant, which we. take as 0.4; d Q [n] is the zero plane displacement height,, 
and Zg[n] is the roughness length. F^[n] is defined as 


F w [n] = 1 - 


9-4Ri B [n] 


l+c|Si B [n] I 1 / 2 ’ 


Rig[n] < 0 (2.7a) 


F w [n] =* 


(l+4.7Hi B [n]) 2 ’ 


0 <Rijj[n] < 0.2 (2.7b) 


where Ri B [n] is the bulk Richardson number and is estimated as, 


Ri [n] = _8-. 2 rM-(TaM-Ts[n]) 

B Wvl 


( 2 . 8 ) 
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where g is the acceleration due to gravity, T a [n] is the air temperature, T s [n] is 
the surface temperature, and z r [n] and V| are the reference height and the 
modified wind speed, and they are expressed following Dickinson et al. (1992) 
and Smith et al. (1993), 

z r [n] = z 2 [n]- d 0 [n] (2.9) 

and 

Va = KC^)) 2 + U c (2.10) 

with U c =1.0 m/s for unstable conditions and U c =0.1 m/s for stable conditions. 
The parameter c in Eq. (2.7a) is expressed as 


c = 49.82 x a2[n] x (^^-^-^ L ) 1 / 2 . (2.11) 

In the Louis (1979) representation, the transfer coefficients for water and heat 
are taken to-be equal, but ..they can be different from the - coefficient for 
momentum which has been shown by experimental and theoretical work (e.g. 
Garratt and Hicks 1973, Garratt 1978, Brutsaert 1982, and Duynkerke 1992). 

Based on the formulation of Blondin.(199l) and Ducoudre et al. (19.93), 
the transpiration rate was obtained using the form 


E t [n] = [!-(• 


Wi[n] 

Wi~R 


2/3 

•) ]E P [n] 


r w| n ] 

r w [n]+r 0 [n]+rc[n ' 


( 2 . 12 ) 


where rc[n] is the canopy resistance. It is expressed as 


r c [n] = 


r mint n l S sm M 
LAI[n,m] 


(2.13) 
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In Eq. (2.13),- r miTl [n] is the minimum canopy resistance, g sm [n] is a soil 
moisture stress factor depending on the water availability in the root zone for the 
nth surface cover class. It is expressed as 


SsL M 

= 1, 


Wj[n]>w[ r 

(2.14a) 

§sm ^ 

' Wj[n] 

‘ w f . 

I 1 

Wj V <Wj[n]<Wj r 

(2.14b) 

®sm W 

- o , 


WjW^Wf 

(2.l4c) 


where Wj[n] is the soil moisture content in layer j, j=l, 2. Wj* is the critical 

value above which transpiration is not affected by the moisture stress in the Soil* 
w 

and Wj is the soil. moisture content at permanent wilting point. Water can be 
extracted from layers 1 and/or 2 depending on fractions of roots f-Jn] and 
f 2 [n] in each layer. 

There is no soil moisture stress, i.e., g sni [n]=i in Eq. (2.14), if either 
(i) W 2 [n] is greater or equal to W 2 , and f 2 [n] > 0.5 or (ii) Wt[n] is greater or 
equal, to W^ , and fjjn] > 0.5. In case (i), the transpiration is supplied by layer 
2 with no soil moisture stress, i.e., E t [fi]=E 2 [n] (regardless of water availability 
in layer 1); • in case (ii), the transpiration takes water from layer 1, i.e., 
EJn]=E|[n], also without any soil moisture stress. Otherwise, the transpiration 
rate is 


E t [n] = fjjn] • E\[h] + f 2 [n] . E£[n] (2.15) 

where E^[n], E|[n] are the transpiration rate from layer 1 and layer 2 

respectively, computed by using Eq. (2.12). If the roots only extend to layer 1, 
then EJn] = Ej[n] with f 2 [n] = 0. 
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For the case of a continuous rainfall with a rate lower than the canopy 
evaporation rate, it is important to consider evaporation from the vegetation 
when there is insufficient intercepted water to supply the atmospheric demand 
within one time step. Thus in general,, the evaporation rate from the canopy 
layer, E c [n), can be expressed as 

E c [n]=f[n]-Ec[n] (2.16) 

where f[n] is the fraction of the time step required for canopy, evaporation to 
exhaust the canopy interception storage. It is given by 


_ Wj[n]+P • At x 

f[n] = min(l, V ; ) 

E c [n] • At 


(2.17) 


where P is the precipitation rate, and At is the time step which is taken as one 
hour in the model calculation. The transpiration during the time step is then 

*w = (10 - -®pW + 


(2 - 18) 


where the first term represents the fraction of the time step for which no 
evaporation occurs from the canopy interception storage, and the second term 
represents the fraction of the time step for which both evaporation from the 
canopy and transpiration occur. 

Evaporation from bare soil is extracted only from layer 1; bare soil 
evaporation E 2 from layer 2 is assumed to be zero. When layer 1 is saturated, it 
evaporates at the potential rate Ep[N+l], i.e., 
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% Ep[N+l] . (2.19) 

When it is unsaturated, it evaporates at rate E x which varies within the bare 
soil area due to the inhomogeneities in infiltration, topography and soil 
characteristics. E r is computed using the Amo evaporation formulation 
(Francmi. and Pacciani 1991). The Amo model uses the structure of the 
Xinanjiang model (Zhao et al. 1980, see also Wood et al. 1992), and assumes 
that the infiltration capacity varies within an area, and can be expressed as 

i/bj 

* - ^m [1 - (1 - A) ] (2.20) 

Where i and i m are the infiltration capacity and maximum infiltration capacity 
respectively, A is the fraction of an area for which the infiltration capacity is 
less than i, and b; is the infiltration shape parameter. . Let A s represent the 
fraction of. the bare soil that is saturated, and i Q represent the corresponding 

point infiltration capacity. Then, as suggested by Fig. 2.2, Ej can be expressed 
as 


El 


As 

= Ep[N+l]{ J dA + j 

0 A s imU-(l-A) ' ‘ ] 


dA}. 


( 2 . 21 } 


la Eq. (2.21), the first integral represents the contribution of the saturated area, 
which evaporates at the potential rate. Since there is no analytical expression for 

the second integral in Eq. (2.21^. Ei is obtained through a power series 
expansion: 


Ei=Ep(N + l|{A 3 + i(l_A s )(l+ bi - • - 1/bi 


J m 

b: , 2/bi b . 3/b: 

■2+E~ (1 - A ») +7^-0 -As) + 


1+^(1 -As) + 

}}. 


( 2 . 22 ) 
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This Approach accounts for the subgrid variability in Soil moisture within the 
area covered by bare soil. 


2.3. Canopy layer 

The water balance in the canopy layer (interception) can be described by 
dWsfnl 

— j^-=P-Ec[n]-P t [a] t 0<Wj[n]<W im tn] (2.23) 

where PJn] is the throughfall rate of precipitation which occurs when Wj^fn] is 
exceeded for the nth surface cover class. 


2.4. Surface runoff from, bare soil 

Surface runoff is computed Using the formulation for infiltration given by 
Eq. (2.20). The Xinanjiang formulation, which is described in detail by Wood 
et id. (1992), is assumed to hold for the upper soil layer only. The maximum 
soil moisture content of layer 1, W l5 is related to i m and bj as follow's, 



*m 

1 + bj * 


(2.24) 


The Xinanjiang model effectively assumes that runoff is generated by those areas 
for which precipitation, when added to soil moisture storage at the end of the 
previous time Step, exceeds the Storage capacity of the. soil. The direct runoff 
from these areas is Q d [N+l], where N+l i ndicates the bare soil class. In 
integrated form, the result is 


Q d [N+l]*At = P-At-Wj +Wj[N+l], i 0 +P-At>im (2.25a) 


QJN+lJ-At = P-At-Wj + Wi[N+l] + Wj [l,iQ t PAt ] 1+b \ 

»m 


i Q +P • At < i m (2.25b) 
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where Wj [N+l] is the soil moisture content in layer 1 at the beginning Of the 
time step. Note that, for the bare soil class, there is no canopy storage, hence 
’’throughfall” is equal to precipitation P. For bane soil, the water balance in 
layer 1 is 


W* [N+l] = WjJN+1] + (P , Q d [N+l] - Q ia [N+l] - Ej) . At (2.26) 

where Wj [N+l] is the soil moisture, content in layer 1 at the end of each time 
step, and Q^N+l] is the drainage from layer 1 to layer 2. Assuming that the 
drainage is driven by gravity, we use the Brooks and Corey (1964) relation to 
estimate the hydraulic conductivity, and thus we can express the drainage rate 
from layer 1 to layer 2 as 


QpIN+l] = K.-t -ffi 1- * ) 


(2.27) 


where Kg is the saturated hydraulic conductivity, 0 r is the residual moisture 
content and can be taken as zero in general due to its very small magnitude, 
and Bp is the pore size distribution index. 


2.5. Subsurface runoff from bare soil 

The formulation of subsurface runoff (baseflow) follows the Arno model 
conceptualization (Francini and Pacciani 1991), which is applied only to the 
lower soil layer (drainage from layer 1- goes only to layer 2, and does not 
contribute to runoff). The baseflow rate is g iven by: 


Q b [N+l] - 


DsDjn 

W S W 2 


■W* 2 [N+1] , 


for 0 < W* 2 [N+1] < W S W 2 (2.28a) 
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Q b t N+1 l = w? w 2 [N+l] + (Dm- 

WgW2 


D S D 


s^m 


W s 


)(■ 


Wr 2 [N+l]-WsW 2 2 
W 2 -W S W 2 * 


for W^[N+1] > W S W 2 (2.28b) 


where Q^[N+1] is the subsurface runoff rate, Djq is the maximum Subsurface 
flow rate, D s is a fraction of D^, W 2 is the maximum soil moisture content of 
layer 2, W s is a fraction of W 2 , with D s < W Si and W 2 [N+1] is the soil 
moisture content at the beginning of the time step in iayer 2. Eqs s (2.28a) and 
(2.28b) describe a recession that is linear below a threshold (Eq. (2.28a)), and 
nonlinear at higher soil moisture values (Eq. (2.28b)) as shown in Fig. 2.3. . The 
nonlinear drainage is required to represent situations where substantial 
subsurface storm flow occurs. Eqs. (2.28a) and (2.28b) have a continuous first 
derivative at the transition from the linear to nonlinear drainage as shown in Fig. 
2.3. 


Using Eqs. (2.28a) and (2.28b), and the notation that W 2 [N+l] is the 
layer 2 Soil moisture content at the end of the current time, the water balance 
for layer 2 is 


W 2 [N+l] = W 2 (N+1] + (Q 12 [N+1] - Q b [N+l] - E 2 ) • M (2.29a) 

when W 2 [N+1] + (Q^IN+l] -Q^[N+1] ~E 2 )*At < W 2 , in which case is 
given by Eq. (2.28a) or (2.28b). _ 

In the case W 2 [N+1] + (Q^N+l] — Q^[N+1] -E 2 ) • At > W 2 (where 
Qb[N+l] is given by Eq. (2.28a) or (2.28b)), 


Wj[N+l] = W§ 


(2.29b). 


and 

Q b [N+l] = W 2 [N+1] + (Q 12 [N+1] - Q b [N+l] - Ej) • At - W 2 . 
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When Eq. (2.29b) applies, the total subsurface runoff from layer 2. is given by 
Q b [N+l] — Q^[N+1] + Q^[N+1]. In practice, this condition occurs rarely. 

2.6. Surface and subsurface runoff from soil with vegetation cover 

The equations for surface and subsurface flow rate, and the water balance 
in each layer are the same for cover classes with Vegetation as for the bare soil 
case, except that P, Ej, and E 2 are changed to P t [n], Ej[n], and E^n], 
respectively in Eqs. (2.25), (2.26), and (2.29), to reflect the vegetation class. 

The total evapotranspiration rate E^ and the total runoff rate Q can be 
then expressed as 

E - £ C v [n] . (E c [n]+E t [n]) + 0 V [N+1] • E, (2.30) 

n=l 

N+l ‘ 

Q = S Cy[n] • (Q Jn] + Q, [n]) (2.31) 

n=l a D 

where C v [n] is the fraction of vegetation cover for the nth (n=l, 2, -• •, N) 
surface cov|p + class of interest, Cy[N+ll the faction of the bare soil covered 
area, and ^ C v [n] = 1. 

nsl 

' 2.7. Aerodynamic flux representation . 

The two-layer hydrological model described above is used in conjunction 
with the energy balance at the land surface and the thermal properties of soils to 
calculate the -surface temperature, and simultaneously, the fluxes of sensible 
heat and ground heat which depend on surface temperature. The energy balance 
equation for an ideal surface of the nth surface cover class can be expressed as 


Rn[n] = H[n] + p w L e E[n] + G[n] 


(2.32) 
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where £[nj is given by 


E[n] * E c [n] + E t [n] , n=l, 2, • • • , N (2^3*) 

E[N+1] * Ej (2.33b) 

where H[n] is the sensible heat flux, /> w L e E[n] is the latent heat flux (e.g., with 
units of Win' 2 ), and G[n] is the ground heat flux. For a surface that is 
relatively flat and homogeneous, the energy balance equation for a layer of the 
air column bounded , by the ground surface at the bottom and a surface of given 
height in the atmosphere above, can be expressed as 


RnW — H[ n ] + P^h^E[n] + G[n] +. AHg[nj 


(2.34) 


where Afig[n] is the change in the. energy storage in the layer per unit time, per 
unit area. The sensible and latent heat fluxes, as well as the net radiation, axe 
associated With, the top surface of the air layer, and the ground heat flux with 
the bottom of the layer. The rate of heat energy 


AH s [nl = 


Ts[n] ) 2 a [n] 
2~At 


(2.35) 


where Tg [n] and T$[n] are the surface temperature of the bottom surface of the 
air layer at the end and at the beginning of a time step respectively, and z a [n] is 
the height of the top surface of the air layer which is used only When AH s [n] is 
considered to be significant. 

The net radiation R n [n] is given by 

ll n [n] = (l-«[ n ])R s + t [n).(R L -<,TsW) (2.36) 

where a(n] is the albedo of the nth surface cover class, Rg is the downward 
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shortwave radiation at the surface, e[n] is the emissivity of the nth surface cover 
class, Rl is the downward long- wave radiation at. the surface, and a is the 
Stefan-Boltzmann constant. The latent heat flux, which is the link between the 
water and energy balances, is obtained from Eq. (2.33). The sensible heat flux 
is given by 


H W = -^-(TsW-.T a [nD (2.37) 

where r^[n] is the aerodynamic resistance for heat transfer. We take r^[n] to be 
equal to r#[n] in Eq. (2.4). The ground heat flux G[n] is estimated through the 
two thermal soil layers (Fig. 2.4). The two thermal soil layers are different from 
the two soil moisture layers (i.e., the upper zone and lower zone) discussed 
above.. For the first soil layer, with soil depth Dj (subsequently assumed to be 
50 mm), we have, 


G W = (T s [n] - Tj[n]) 


(2.38) 


where «[n] is the soil thermal conductivity, and TjJn] is the soil temperature at 
depth Dj. For the second .soil layer with depth D 2 , at which, the bottom 
boundary condition is a constant Soil, temperature T 2 , the law of energy 
conservation (assuming that the heat storage . in the first soil thermal layer is 
negligible) gives, 


a(C s [n)T[nj) _ eH f [n] 

at w~ 


(2.39) 


where C s [n] is the soil heat capacity, T[n] is the soil temperature, H^nj is the 
heat flux, and t and z are the time and the soil depth respectively. Assuming 
that C s [n] does not change with time, from Fig. 2.4, Eq. (2.39) can be written 
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as 


CsW-^|M-= __ G”M-G’[a] 


(2.40) 


where AT[n] is the .temperature gradient, G”[n] is the heat flux across the 
bottom boundary at depth D 2 , and the G [n] .is the heat flux across the soil face 
at depth (see Fig. 2.4). Since it is assumed that there is no heat Storage 
within soil depth Dj, we have 


G[n] = G[n] . 


(2.41) 


. 4 11 

In addition, G [n] can be expressed as 

g ”M = -^( t iM-t 2 ) 


(2.42) 


If we assume that 


T[n] = Tlln l +T2 


(2.43) 


then 


AT[n] _ TiM-Tj[n] 
At 2 m 


(2.44) 




where [n] and T^[n] are the soil temperature at depth Dj at the end and the 
beginning of a time step respectively. 

By substituting Eqs. (2.41-2.44) into Eq. (2.40), we can obtain 


CsfcHT! [n]- Tj[n) ) G[n] 
2 • At D 2 


- T2) (2.45) 

d 2 


At present, C s [n] and «[n] are not considered to be functions of the soil water 
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content (although such an adjustment would be straightforward), and are taken 
to be the same for both soil thermal layers. From Eqs. (2.38) and (2.45), the 
ground heat flux G[n] can be expressed as, 


G[n] = 


^2 


(TsN-T 2 ) 


1 + 


Di 

~T57 


+ -g ?M^2 (T.M-ntnl) 

. C s [nj • Dj • P 2 ~ 

2 • At • ri[n] 


(2.46) 


For the case where AH s [n] is not significant or there are not enough data 
available to evaluate the energy balance within a layer, the energy balance 
equation for an ideal surface (Eq. (2.32)) can be used instead of Eq. (2.34). 
From Eqs. (2.36), (2.37), (2.46), and Eq. (2.33) (scaled by the latent heat of 
vaporization and the density of liquid water), we can obtain the sensible heat 
and ground heat fluxes and the surface temperature, for the nth cover class. In 
the case where AH s [n] is negligible, the surface temperature T s [n] is solved 
iteratively from Eq. (2.47) below, 


e[n](xT s [n] + (• 


^a c P 


+ 


T57 + 


Cg[n] 'D 2 
2 At 


1+ 


D 


1 


+ 


n 


• Dj -D 2 

2 * At • k[u] 


) • Ts[n] = (1 - a[n] )Rg + 


[n]R L + T a [n] -,„LeE[n] + 


K ( n ] ‘^2 , Cg[n] -D 2 -T^[n] 
Do 2 * At 


1+ T57 + 


D i • D2 


(2.47) 


2 • At • «[n] 


For the case where AH s [n] cannot be ignored, Eqs. (2.34) to (2.37), and (2.46) 
axe combined to give 
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"■ + r„tt4 / p a c P ^a c P 2 aM ^ ^2 

I 

1 + 


*M ' CsM • p 2 


+ 


.WKTsJnir + + T~ bl ' . ' Csi^-D, ^ W - 


T5T + 


2 - At 
2 • At • /c[n 




+ <W-&l+ Ag-T a W- fw L e EH + . W ^ t T,H 


+ 


/t[n]-T 2 C s [n]-D2-tj[n] 
Do 


2- At 


D i. 

h 


1+ “ST + 


Cg n • D ^ • D 2 
2 • At • *[n] 


(2.48) 


T s [n] is determined in the same manner as for Eq. (2.47). Therefore, the 
effective surface temperature T s , sensible heat flux H, and ground heat flux G 
can be obtained as, 


N+l 


7! C v [n] -T s [n] , 
n=l 

(2.49) 

N+l 

£ C v [n]-H[n] , 

n=l 

(2.50) 

N+l 

£ Cv[a].G[n] . 
n=l 

(2.51) 


2.8. Snow 

When snow is present, the model is. coupled with a single-layer, energy- 
and mass-balance snow accumulation andablation model (Wigmosta et al. 1994). 
At the snow-air interface, the.energy exchange is described by the net radiation, 
sensible heat, evaporation from the water in the snowpack and sublimation or 
condensation, and the heat, advected to the snowpack by rainfall. The snow- 
ground interface is assumed to be a zero ener gy flux boundary . Snow albedo is 
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determined based on snow age. The present version of the snowmelt model does 
not consider fractional snow coverage; it is assumed that the entire area is 
covered, by a uniform depth of snow if a snowpack is present. 

2.9. Model calculation procedure 

In the formulation of the two-layer VIC model, many variables are a 
function of the surface temperature. For example, the surface temperature is 
needed to calculate the bulk Richardson number, vapor pressure deficit, and 
net radiation. Once the bulk Richardson number is determined, the 
aerodynamic resistance can be calculated, accounting for stability correction. 
Given the stability-corrected aerodynamic resistance, the Vapor pressure deficit, 
and net radiation, the potential evaporation can then be estimated. However, 
computation of the surface temperature, requires an iterative solution of Eq. 
(2.47) or (2.48), which is implemented as follows: 

(1) . Set the surface temperature to the air temperature .at the first time step. 
This allows computation, of the initial values of the bulk Richardson number, 
vapor pressure deficit and net radiation that are needed to estimate Ep[n] 
through the Penman-Monteith formulation. 

(2) . Iterate Eq. (2.47) or (2.48) to solve for the surface temperature. 

(3) . Use the surface temperature obtained from step (2) to calculate the bulk 
Richardson number, vapor pressure deficit, and net radiation again. 

(4) . Recalculate the surface temperature iteratively using Eq. (2.47) or (2.48). 
The surface temperature obtained from this step is then considered to be the 
surface temperature of the first time step of the model simulation. 

(5) . -For. subsequent time steps, use the surface temperature- from the previous 
time step to calculate the bulk Richardson number, vapor pressure deficit, and 
net radiation, then repeat steps 2-4. 
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The procedure described above is not iterative in the same sense as the 
procedure used to solve for the surface temperature from Eq. (2.4?) or (2.48), 
since the steps are only repeated once. The use of a single -iteration is justified 
by the relatively smooth variation usually observed in surface temperatures due 
to the thermal inertia of the soil column. Of course, multiple iterations could be 
performed if required. Such an approach in fact implies two nested iterations; 
one to solve Eq. (2.47) or (2.48), and the other to determine the bulk 
Richardson number and related quantities needed to compute the surface energy 
fluxes. 

2.10. Summary 

A generalization of the VIC (Variable Infiltration Capacity) model which 
incorporates a twO-layer description of the soil . column has been described. . In 
the' soil. Column, the upper layer is characterized by the usual VIC spatial 
distribution of soil moisture capacities, and the lower layer is spatially lumped 
and uses the. Amo (Francini and Pacciani 1991) drainage representation. The 
model is designed for application Within coupled land-atmosphere-ocean GCMs, 
such, as are used for numerical weather, prediction and global climate simulation. 
The model partitions the area of interest (e.g., grid cell) into N+l land surface 
cover types; for each land cover type the fraction of roots in the upper and lower 
zone is specified. Evaporation occurs via canopy evaporation, evaporation from 
bare soils (land cover class-N+1) and transpiration, which is represented with a 
canopy and architectural resistance formulation. 

The two-layer VIC model described here has been tested against observed 
data. The behavior of the model is described in Chapter 3 using surface flux 
data from two sites: the FIFE experiment in central Kansas, and the 

ABRACOS experiment in Brazil. 
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CHAPTER 3 MODEL APPLICATION 


In Chapter 2, the two-layer VlC model, was described. In this chapter, 
the performance of the model is evaluated for two applications. 

3.1. FIFE site 

The first test location for the model was the FIFE (First ISLSGP Field 
Experiment) site in central Kansas in the United States. The FIFE site is a 
15x15 km2 region on the Konsa Prairie, a native grassland preserve near 
Manhattan, Kansas.. It has a fairly homogeneous tall grass cover. The Kings 
Creek catchment, of area .11.7 km 2 , lies within the FIFE site. The FIFE site is 
of interest because of the detailed, measurements of surface fluxes that were 
collected in the summer of 1987. A detailed description of the site is provided by 

Sellers et al. (1992). 

3.1.1. Data .description 

During the period MayOctober, 1987, four intensive field campaign 
(IFCs) were conducted at the FIFE rite, during which tower-baaed 

measurements of latent, sensible, and ground heat fluxes were made (Sdlers et 

al. 1992). In addition, throughout the summer of 1987, a network of porta e 
automated meionet (PAM) stations was.operated.. from which meaaured values 
of i r-— i*g solar and long-wave radiation, and other meteorological data are 
available. Furthermore, long-term streamflow data exist for Kings Creek, along 
with long-term dimatalogical data at nearby Manhattan, KS, which allows for 
validation of the hydrologicri portion of the VIC model. The strategy for 
validation of the model was to estimate the hydrological parameters using 
precipitation and Streamflow data for part of the long-term Kings Creek ^record, 
Ld to evaluate its hydrological performance using the remaining part of e 
record. The model’s surface flu* algorithms were then parameterized, an 
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validated using measured fluxes observed during the summer, 1987, IFCs. A 
schematic of the FIFE site with approximate locations of King’s Creek 
catchment (shaded area), flux stations, and meteorological stations is shown m 

Fig. 3.1. 

Daily precipitation and temperature maxima/minima have been collected 
at Manhattan, KS, which is about 11 km from the centroid of the Kings Cree 

-set, . since the late 1800s. Daily average stream discharge data for Kings 

Creek (U,S. Geological Survey Station No. 06879650, 11.7 km 2 ) have been 

collected since about 1980. Surface meteorological and surface flux data at the 
FIFE site are limited to selected periods during the summer of 1987. Data from _ 
the PAM stations include surface pressure (p), mimng . ratio (w) and air. 
temperature (T») at the 2 m level and horizontal wind speed (u) measured 5.4 m 
above ground level, surface temperature (T.). ground sod temperature, Tu 
and Ten, at 10 cm and 50 cm below the surface, respectively, and. downward 

^ j. . p j- ai-ri data were collected from flux 

short and long-wave radiation. Radiation data were 

stations (eddy correlation and Bowen ratio). Data from both PAM **»•'£ 
flux stations were averaged for each date and time among all the s *»* j* 
Betts et al. (1993). They found from consistency analysis of the calculated 
measured net radiation that the flux data were more self-consistent than the 

PAM data. Therefore, we used the radiation data from the flux stations and the 

atmospheric data 60 m PAM stations to test our model surface flux an . ace 
temperature predictions. Data for 35 days common.to the two data sets m th 
summer of 1987 were used. They are June 30 -July 11, August 9-20, and 

October 6-16. 


3.1.2. Parameter estimation 

The model parameters can be classified into hydrological parameters and 
atmospherically related parameters. The hydrological parameters include the 
infiltration shape parameter bj (Eq. (2.20)), the mil pore sire distribution nde 
Bp the residual moisture content e t , the saturated hydraulic conductmty Ks 
(Eq. (2.27)), the three base flow-related parameters D„ D a , and (Eq- 
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C 0 

(2.28)), the maximum soil moisture contents Wj 'and W 2 iii layers 1 and 2, 
respectively (Eqs. (2.24) and (2.28)), the wilting point wj* and the critical point 
Wj (j=l, 2) in Eq. (2.14). Atmospherically related parameters include 
architectural resistance rg[n] (Eq. (2.1)), minimum canopy resistance r^^n] 
(Eq. (2.13)), leaf area index LAI[n,mj (m=l, 2* • • •, N; m=l, 2 4 . • 12) for 

each surface cover class (Eq. (2.13)), the Zero plane displacement height dg[n], 
roughness length Zq[u], and the relative fraction of roots in each of the two soil 
zones fjn] and f 2 [n] (Eq. (2.15)). We classify fy[n] and f 2 [n] as atmospherically 
related parameters because they determine the canopy resistance (Eqs. (2.13- 
2.15)). 

Among the hydrological parameters, only three (bj, D s , and W s ) are 
best estimated using streamflow data if they are available. (both.A s and ig in 
Eqs. 2.22 and 2.25 are not model parameters, they are evaluated at each, time 
step).. The other hydrological parameters can be estimated using soil 
characteristics. Clearly, for application in GCMs, global, parameter estimation . 
using streamflow data is infeasible; for GCM applications. Dumenil and Todini 
(1992) have suggested values for bj, Dg, and .Ws- An ongoing research topic, 
which will be investigated in the GEWEX . Continental Scale International 
Project (GCIP) is to develop regional relationships for GCM. hydrological 
parameters. However, because streamflow data were available for Kings Creek, 
we made use of the observed data to estimate bj, Dg, and W s . 

Since the FIFE .site has a fairly homogeneous, tall . grass cover, the 
number of vegetation types n in this model application can be taken to be 1, 

1. e., C v [n]=:1.0. To estimate bj, D s , and Wg through calibration, we need to 
know Ep[l] md r w [l], in addition to parameters Wf, Wj 7 , Wj r (j=l, 2), K s , 
B P> *r, fy[l]> f2[l]. u n (z 2 ), d 0 [l], z 0 [l], r 0 [l], r mi - n [l], and LAI[l,m] (m=l, 

2, 12). At the FIFE. site, the data required to estimate Ep[l] by the 
P enman-Monteith method and r w [l] by Eq. (2.4) are available only during the 
IFCs. Therefore, for the purposes of estimating the hydrological parameters, 
we used Hamon’s method (Hamon et al. 1954, Hamon 1961) which requires only 
daily air temperature and latitude to estimate Ep[l]. During summer, 1987 we 
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also used Hamon Ep[i] during the between-IFC periods to allow continuous 
computation of soil moisture (needed as initial Values during the IFCs). The 
daily Ep[l] computed using the Hamon formula was compared with the daily 
Ep[l] obtained using Penman-Monteith’s equation for the 35 days of the 1987 
FIFE IFCs. The comparison indicated that the Hamon equation gives smaller 
Ep[l] estimates, but the pattern over the 35 day period was s imil ar for both 
Ep[l] estimates. Therefore, we scaled the Hamon estimates to have the same 
mean as the Penman- Monteith estimates, using an adjustment factor kg, which 
was determined to have a mean of 1.64 with a standard deviation of 0.70. The 
scaled Hamon estimates were used for the long-term hydrologic water balance 
computations, except during the IFC periods, when the data needed for 
computation of the PenmamMonteith Ep[l] were available. During the IFC 
periods, Ep[l] and r w [l] Were estimated by the Penman-Monteith method and 
by Eq. (2.4), respectively. 

For the 35 days of the IFCs, we calculated an average aerodynamic 
resistance (equal to the inverse of the. product of the drag coefficient from Eq. 
(2.6) and the wind speed under the assumption that the resistance to the transfer., 
of momentum, and .water are equal). This average aerodynamic resistance was 
then used for the purpose of estimating the hydrological model parameters, and 
for computing the soil moisture at the be ginning of the first. IFC and between 
IFCs (but not for validation of the energy fluxes during the IFCs reported. in 
Section 3.2). The average aerodynamic resistance over the 35 days was 40.8 s/m 
with a standard deviation of 29.7 s/m. This value is within the range given for 
short grass and crops by Monteith and Unsworth (1990). Since the roughness 
length of many crops decreases as wind speed increases, the inverse of 
aerodynamic resistance is approximately a constant over a range of .low wind 
speeds. The daily average wind speed during the 35 days was 2.38 m/s, and the 
aerodynamic resistance (40.8 s/m) was taken as constant for the estimation of 
the three hydrological parameters. In addition, we did not correct for 
atmospheric stability, primarily to assure compatibility of Ep[l] between the 
IFCs and during the longer period of hydrological water balance simulation, 
when the data needed to make the corrections were not available. However , the 
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stability correction given by Eqs. (2.4) and (2.5) was applied to the energy flux 
computations performed for the model validations reported in Section 3.2. 

At the FIFE site, the depths of layer 1 (upper zone) and layer 2 (lower 
zone) Me about 0.5 m and 2.5 m respectively (Famiglietti and Wood 1993). 
Since the soil texture at FIFE is silt loam (EPA 1991), the porosity was taken 
to be 0.5, and thus W^=0.25 m and W2=1.25 m* and Wj 7 and Wj" r are about 
26% and 46%, respectively, of the total water that the soil can hold. However, 
Smith et al. (1993) reported that evapotranspiration was not observed to be 
limited by soil moisture in the 20% -30% range, and. they took 18% as the 
wilting point instead, which we also used as our estimate. .In this study, we 
used 70% of field capacity as our critical point (we found via sensitivity analysis 
that , almost the same results were obtained when .the critical point Was 75% . of 
field capacity). The K s , Bp, 0 r were- taken as 6.44 mm/h, 0.16, and 0.01 m 
respectively, following Famiglietti and Wood (1993). Since the vegetation is 
dominated by grass, we assumed that all the roots are .in the upper zone (i.e., 
^[11=1.0 and f 2 [l]=0.0). 

Because the wind speed from the PAM stations Was measured at. 5.4 m 
above the ground surface, and the other meteorological data were measured , at 
Z 2 [l ]=2 m above ground surface, the wind speed was converted to the 2 m level 
through a logarithmic velocity profile. Sugita and Brutsaert (1990) estimated 
the zero plane, displacement height do(l]=26.9 m, and the surface toughness 
length Zq[1]=1.05 m at FIFE by analyzing neutral . wind velocity profiles 
measured by radiosondes. They found that a logarithmic velocity profile only 
holds over the height ranges between 50 m ± 19 m and 202 m ± 101 m above the 
ground surface. However, their values should be interpreted in the context of 
the Flint Hills region, which is characterized by relief of about 25m between 
steep ridges and valleys. By contrast, Smith et al. (1992a) used much smaller 
local values of d()[l]=0.25 m, and Zq[1]=0.07 m. Their values fall between uncut 
grass and long grass/crops for a relatively flat area (Ary a 1988). Since the FIFE 
site is only a small part of the Flint Hills region which covers a 50-80 km wide 
north-south strip in Kansas from Nebraska to Oklahoma, we decided to use the 
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smaller values for dg[l] and Zq[1], and assumed a logarithmic velocity profile 
locally. The 2 m wind speed can then be estimated as 


u n (z 2 )=u n( , 1 )-— lf 0' d( (3.1) 

M ~ olH ' > 

where. z^[l]=5.4 m, and u n (z^) is the corresponding measured wind speed (m/s). 
The value of r Q [l] for grassland is taken as 2.0 s/m (Ducoudre et al. 1993). 
Monteith and Unsworth (1990) suggest that for crops r jnin [l]=l00 s/m. Smith 
et al. (1993) found, based on optimization, that the best values of r min [l] are in 
the range 100-125 s/m. We take r fnia [l]=100 s/m to be consistent with 
Monteith and Unsworth (1990)._ 

• The monthly average LAI[l,m] (m=*l, 2, • • •, 12) were derived from the . 
average normalized difference vegetation index . (ND Vis) given by EPA (1991) . 
with LAlmax=6.0 and LA I m ^ =0.1 which are consistent with the values used by 
Smith et al. (1993), 

LAI[l,m] = 0.1 + 0.0628(NDVI[l,m] - 53.0). (3.2) 

The average monthly NDVIs for 1986, 1987, and 1988 at FIFE are listed in 
Table 3.1. 

The hydrological parameter D m can be either estimated by identifying 
extended dry periods during the calibration interval 1982-85 - using the 
precipitation data, and recession rates inferred from the observed Kings Creek 
streamflows during, these periods; or by multiplying saturated hydraulic 
conductivity by an average soil slope. We used the first approach, which gave 
Dm = 8.2 mm/d. The hydrological parameters bj, D s , and W s , were 
estimated using streamflow at Kings Creek, and precipitation, and 
maximum/minimum temperature data at Manhattan, KS from 1982-1985. The 
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calibration gave b^O.008, D s =7.7xlG' 5 , and W s =0.96. The one-layer 

snowmelt . model was not used to obtain the above model parameters, since not 
much Snow occurs in the Kings Creek catchment. Hydrographs for two of the 
calibration years (1983 and 1984) are shown in Figs. 3 . 2 a and 3.2b. The model 
reproduces the streamflow reasonably well; discrepancies are attributed to a) 
the distance of the precipitation gage from the Kings Creek catchment; b) the 
inability of a single gage to represent spatial variations in precipitation; c) the 
use of a daily time step for a. relatively small catchment whose time of 
concentration is on the order of an hour or less; and d) small scale 
heterogeneities which can strongly affect runoff production in small cat chm ents, 
and Eire not captured by a macroscale model such as VIC. 

With the parameters described above, together with the parameters a[l], 
*[!]> C s [l], D 2 , and Tg, we then used the PAM and flux data to test our 
model-predicted surface fluxes and Surface temperature against the measured 
values. The albedo a[l] was taken as. 0.2 during the IFCs following Famiglietti 
and Wood (1993). The thermal conductivity «[ 1 ] and soil heat capacity C$[l] in 
Eq. (2.46) were estimated to be 0.514 Wm-^-K * 1 and 2.13 x 10 ** Jm-^K' 1 , 
respectively, following Smith et al. (1992b, 1993). The depth D 2 was taken, to 
be 0.45 m, and the temperature T 2 (i.e., Tgg) in Eq. (2.46) was prescribed as 
293.6 °K, which was the average of T 5 Q for the selected 35 days of the IFCs. 
The standard deviation of T 5 Q for the 35 days was 3.1 °K. The values of the 
hydrologically and atmospherically related model parameters are listed in Table 
3.2. We compared the surface energy budgets computed using both Eq. (2.32) 
and Eq. (2.34), and found that there was almost no difference in the results 
when we took 3 a [l]=Z 2 [l ]=2 m. 

3.2. Model validation at FIFE site 

Figs. 3.3a and 3.3b show predicted and observed streamflow for 1986 and 
1987, two years not in the calibration period. Generally, the results are 
consistent with those of the calibration period — the dry period flows are fairly 
well represented, as is_the timing of the major peaks, but the magnitudes of the 
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peaks, especially the largest ones, are subject to major errors. . In this study, 
streamflow prediction is not of primary importance; the purpose in evaluation of 
the predicted hydrographs is to provide evidence that the model is producing a 
plausible soil water balance. To this extent, the hydrograph simulations were 
adequate. 

After, estimating the hydrological model parameters, we used the FIFE 
surface fluxes and meteorological measurements for. the summer of. 1987 to test 
the model predictions of latent heat, sensible heat and ground heat fluxes, and 
the surface temperature. We used the Kings Creek, precipitation network, as 
Well as the precipitation, air temperature, and downward solar and long-wave 
radiation composited from the PAM and flux stations by. Betts, et al. (1993) to 
test the model heat fluxes and surface temperature. Results are shown in Figs. 
3.4, 3.5, and 3.6 for parts of the June, July, August, and October IFCs. 

Fig. 3.4 shows part of IFC 2 (from June 30 -July 11). There were 
precipitation events on June 30 and July 7. - On the rest of the days, there was 
little or no rainfall. During this period, the latent heat -flux for dry days was 
typically about 400 Wm~2. The model predicted the latent heat and sensible 
heat fluxes fairly Well, except that it somewhat underpredicted the July 9, 10, 
and 11 latent heat fluxes and overpredicted the sensible heat fluxes on the same 
days. These days were characterized by relatively high winds, high potential 
evaporation, and high soil moisture. The modeled surface temperatures agree 
with the observed ones quite well, but the magnitude of the diurnal cycle of the 
ground heat flux was underpredicted on some of the days. 

Fig. 3.5 shows predicted and observed latent, sensible, and ground heat 
fluxes, and surface temperature, for the August 9-20, 1987 portion of the third 
IFC. Rainfall occurred on August 12, 13 and 18. Before the August 12-13 
storm, the soil was moderately dry. During this period, the observed latent 
heat fluxes were less than 300 Wm*2. After the rainfall, the latent heat fluxes 
increased to about 400 Wm*2,_ During this period, the model predicted the 
latent heat and sensible heat fluxes quite well, except during the nights of 
August 14 and 15, when the latent heat fluxes were overpredicted and the 
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sensible heat fluxes were underpredicted. This is mainly due to the high 
potential evaporation obtained during that time. From Eqs. (2.1) and (2.12), it 
can be seen that large evaporation Would be obtained if the potential evaporation 
is large, even though r w [n] and r c [n] are reasonable. During this period, the 
ground heat fluxes were predicted reasonably well, although the magnitude of 
the diurnal cycle Was underestimated. The surface temperatures were well 
predicted in general* 

Fig. 3.6. shows the energy fluxes and surface temperature for October 6-16, 
1987, a portion of IFC 4 which was characterized by low soil moisture. During 
this period, the observed latent heat fluxes were about 100 Win' 2 or less, while 
the sensible heat fluxes increased to about. 300 Wm -2 (from about 200 Wm' 2 in 
July and August). The model predicts the latent and sensible heat fluxes, and 
surface temperature reasonably well, but it overpredicts the ground heat fluxes 
on most of the days during this period. 

In general, the. model performed satisfactorily, - especially, given its 
simplicity. There are some caveats in interpretation of the results. First, the 
FIFE site is a native grassland, which is characterized by a single vegetation 
type. Therefore, the portion of the model dealing with heterogeneous vegetation 
was not exercised in these tests, so . the effects of certain associated 
simplifications are not reflected in the results. A second, related limitation is 
that since the FIFE vegetation is all grassland, the algorithms dealing with 
trees, which usually extract moisture from the lower, rather than the upper, 
soil moisture zone have not been exercised. The model has, however, been 
implemented for a tropical forest application in connection with the Project for 
Intercomparison of Land-surface Parameterization Schemes (PILPS) (Pitman et 
al. 1993, Liang et al. 1993) and the model results were comparable to those of 
most of the participating models... 

3.3. ABRACOS site 

The second application was to the field site of the Anglo-Brazilian 
Amazonian Climate Ob c ervation Study (ABRACOS). The ABRACOS site is 
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located in a ranch clearing surrounded by forest in the moist, tropical rain forest 
climate regime typical of Amazonia. The site was selected to characterise the 
local and large-scale effects of the amazonian deforestation that has taken place 
over the past 20 years (Myers 1091). ABRACOS is an ongoing project Which 
started in 1990 and will continue through 1994 (Shuttleworth et al. 1991). The 
site is extensively instrumented for micrometeorological and climatological 
measurements. During the three year period of routine monitoring of the near- 
surface climate and soil moisture, five intensive micrometeOrological and plant 
physiological studies are planned to be carried to supplement the long-term . 
routine measurements. 

The time frame of the intensive measurement campaigns is summarized . 
in Fig. 3.7 (taken from Shuttleworth et al. (1991)). One of the objectives of 
ABRACOS is to collect data from cleared Amazonian forest and thus to provide 
energy and. water balance and near-surface climate measurements- for GCM 
studies. The clearing site selected in the ABRACOS project is at. Fazenda 
Dimona, which is located about 100 km north Manaus. (2°19’S, 60°19’W) in 
central Amazonia (Fig. 3.S). It is a typicaLlarge cattle ranch created by felling 
and burning the primary . forest and sowing the clay soil With hardy pasture 
grasses about 12 years ago. The studied area consists of about 84% grass, 11% 
bare soil, 5% trunks, and less than 1% bushes. The height of the grass cover 
Was about 28 cm in September 1990. This selected Site was well managed, and 
there was no overgrazing. A more detailed description of the site is given by 
Wright et al. (1992). The climate is moist and hot with mean a.fltui a.1 rainfall of 
about 2400 mm. The driest months jure from July to October, and the wettest 
months are from March to May. 

3.3.1. Data description 

During 1990 and 1991, two. intensive field measurements, Mission 1 and 
Mission 2, were carried out. Mission 1 started on September 15 and ended on 
November 5 in 1990, While-Mission 2 was from June 30 to September 11 in 1991. 

Four major instruments were used for micrometeorological data collection 
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at . the site. These are a profile tower with six logarithmically spaced 
anemometers and psychrometers, a Hydra eddy-correlation device, a Bowen 
ratio device (from Campbell Scientific Ltd.), and an automatic weather station. 
Net all-wave radiation, specific humidity, and horizontal wind speed were 
recorded at the profile tower. Soil temperatures were recorded by thermistors 
located 15 m .upwind of the profile tower. Ground heat flux was measured as the 
average of nine soil heat flux plates that were installed at 5 mm depth under the 
grass and bare soil, located at equal distances along a transect about 15 m 
upwind of the profile tower. The Hydra device, sampled with a frequency of 10 
Hz, recorded air temperature, specific humidity, and wind speed. A net all- 
wave radiometer was connected to the Hydra. The Bowen ratio device measured 
the temperature and humidity gradients between levels of 0.9. m and 3.2. m. In 
addition, it had its own net all- wave radiometer. The automatic weather station 
recorded the wet bulb . temperature, dry bulb temperature, net all-wave 
radiation, wind speed, total incoming shortwave radiation, reflected shortwave 
radiation, soil heat flux, and rainfall. Hourly rainfall was measured by a. 0.2 
mm tipping-bucket rain-gauge attached to the automatic weather station 
(Wright et al. 1992). All of the four, instruments were located within a few 
meters of each other (Shuttleworth 1993, personal communication). 

The hourly summary data for the first iwo ABRACOS missions Were 
provided, by INPE (Instituto Nacional de Pesquisas), Sao Paulo, Brazil. These 
data Were aggregated from ten minute measurements and Were subjected to a 
data quality control process by staff at INPE. For purposes of determining 
consistency between measurements , the following priorities were, assigned: 1) 
profile tower, 2) Hydra eddy-correlation device, 3) Bowen ratio device, and 4) 
automatic weather station. Missing and unreliable data were indicated by 
” — 99” in the data base. The two-layer VIC model validation was conducted 
using 17 days from Mission 1 and 42 days from Mission 2 that hacLcontinuous 
data. 

Mission 1 started on September 15 and .ended on November 5, but there 
were only 18 days (from October 15 to November 1) that have almost continuous 
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measurements of net radiation, precipitation, Specific humidity, air 
temperature, and wind speed that are needed to drive the two-layer VIC model, 
and latent heat flux, sensible heat flux, and the ground heat flux that are 
needed to validate the model predictions. Soil temperatures Were measured- at 
depths of 0.05 m, 0.1 m, 0.2 m, and 0.4 m. Soil moisture was measured weekly 
in general, and twice each week during the two missions using neutron probes. 
A schematic outlining the profile of the soil water measurements is shown in Fig. 
3.9. The soil moisture in the soil surrounding the eight neutron probe tubes was 
measured at the depths of 0.1 m, 0.2 m, 0.4 m, 0.6 m, 0.8 m, 1.0 m, 1.2 m, 
1.4 m, 1.6 m, 1.8 m, and 2.0 m below the ground surface. Since soil moisture 
was not measured on October 15, but on October 16, only 17 days (from 
October 16 to November 1) were finally selected. Among the selected 17 days, 
the wind speed at 0200 LMT on October 16. (Julian day 289) was missing, and 
was estimated by interpolation. In addition, a few -wind speed measurements on 
October 16 were taken from, the automatic Weather station at the elevation of 2.0 
m above the ground. These wind speeds were corrected, to the 3.6 in level so 
that they were consistent with the wind speed in. the rest of the data ox the 
selected 17 days. The individual flux data Were checked for errors to make sure 
they fell within the ranges given by Betts et al. (1993) as follows: 


Solar radiation 
Solar reflected 
Net radiation 
Latent heat 
Sensible heat 
Ground heat 


-5 to 1200 (Wm* 2 ) 
-5 to 250 (Wm* 2 ) 
-98. to 1000 (Wm* 2 ) 
- 100 to 500 (Wm* 2 ) 
-200 to 500 (Wm* 2 ) 
-150 to 300 (Wm* 2 ). 


The check against the above criteria indicated no gross discrepancies in the flux 
data. -However, when the measured net radiation was compared with the sum of 
the derived latent and sensible heat fluxes, and the measured ground heat flux, 
three hourly values were found to have large relative errors. These are LMT 
2100 on Julian day 293, 1900 on Julian day 301, and 0700 on Julian day 303 
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with relative errors of 176%, 198%, and 41% respectively. These data Were 
replaced by interpolated values. The corrected data agreed well with the 
consistency check. The average wind speed for the 17 days Was 1.67 m/s with a 
standard deviation of 1.10 m/s; the average of the specific humidity was 17.5 
g/kg with a standard deviation of 1.0 g/kg; the average of air temperature Was 
26.6°0 with a standard deviation of 3.7°C. The total precipitation of the 17 
days Was 6.6 mm, indicating that the period Was generally dry. 

The surface air pressure Was not measured during either Mission 1 or 
Mission 2. However, there Were 26 days (from April 18 to May 13) in 1987 that 
had surface pressure measurements at GMT 0000 on each day, and 5 days (from 
May 1 to May 5) in 1987 that had six pressure measurements on each day at 
QOOO, 0600, 1200,. 1500, 1800, and 2100 GMT. The first , set of pressure 
measurements (26 days) had an average of 1004.1 mb with a standard deviation 
of 1.1 mb. The second Set of measurements (from May. 1 to May 5) had an 
average of 1003.5 mb with a standard deviation of 1.6 mb. Based on these 
limited available pressure measurements and discussions with INPE scientists, 
for the model. runs, the surface pressure was taken to be constant 1003.8 mb, 
which is an average of all the measurements. 

Mission 2 Started on June 30 and ended on September 11, 1991, and 
there were some missing observations. Among the 74 days, 72 days (from June 
30 to September 9) have almost continuous measurements of the quantities that 
are required to drive the two-layer VIC model and to compare the model 
predictions .with the observations. The Betts et al. validation criteria for . solar 
radiation, solar reflected radiation, net radiation, latent heat, sensible heat, 
and ground heat fluxes were applied as they Were for the Mission 1 data set. The 
results indicated that all the data were within the specified ranges. However, 
the data passed the consistency check for net radiation on only 42 days (from 
July .5 to August 15). Large errors occurred on other days as shown in Fig. 3.10. 
These large inconsistencies between the observed net radiation and the calculated 
net radiation from the derived latent and sensible heat fluxes, and the measured 
ground heat flux at the end of Mission 2 (i.e., after August 15) were probably 
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due to drift in the hygrometer (Wright 1993* personal communication) . 
Therefore, only the 42 days (from July 5 to August 15) in Mission 2 were Used 
for model testing. 

Among the 42 selected days, there Were six measurement s of specific 
humidities at 1.2 m height, while the rest were measured at 3.6 m. Since the 
specific humidities do not change much between the heights of 1.2 m and. 3.6 m 
above the surface, the 1.2 m measurements were not corrected for height. The 
average Wind speed for the 42 days was 1.43 m/s with a standard deviation of 
0.98 m/s; the average specific humidity was 16.5 g/kg with a standard deviation 
of 1.1 g/kg; the average air temperature was 24.8°C with a standard deviation of 
3.3°C; and the total precipitation of the 42 days Was 105 mm. Although still 
relatively dry in a cliinatalogical sense, the Mission 2 data include several storm 
periods. From Fig. 3.14 it can be seen that there were major storms before July 
5. 

3.3.2. Parameter estimation. 

The ABRACOS site consists of about 84% short grass with average height 
0.28 m, 11% bare soil, 5% tree, trunks, and less than .1% bushes (Wright et al. 
1992). Since measurements were not conducted under the surface covers of the 
trunks-and bushes, these trunks and bushes are assumed to be bare soil which 
then covers 16% of the -area. This consideration is consistent with the 
measurements of the soil heat flux described by Wright et al. (1992). Thus, the 
number of Surface covers at the clearing site is two, with C v [l]=0.84 (n=l for 
grass) and C v [2]=0.16 (N+l— 2 for bare soil). 

As described in Chapter 2, the two-layer VIC model has two soil layers, 
an upper zone and lower zone. The upper zone in the model was designed to 
represent the dynamic behaviour of the soil responding to rainfall events. In 
other words, the soil moisture in the upper zone varies dynamically with rainfall 
events and atmospheric moisture transport conditions, while the soil moisture in 
the lower zone tends to characterize the seasonal soil moisture variations. The 
lower zone only responds to rainfall when the upper zone is relatively wet and 
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thus it can separate the subsurface flow from the rainfall quick response. At the 
clearing site, the soil moisture measurements at the eleven different soil depths 
described above showed that the soil moisture storage in the first 1 m depth 
varied between 0.3 m and 0.45 m from mid-September 1990 to December 1991, 
while th moisture storage between 1 m and 2 m depth varied only from about 
0.45 m to about 0.5 m for the same period (Hodnett et al. 1993). Therefore, the 
upper layer was taken to be from the surface to a depth of 1 m. Following 
Hodnett et al. (1993), the depth of the lower layer was taken to be 9 m. 

The soil at the clearing site consists of clayey oxisols from the Tertiary 
Barreiras sediments. According to the Brazilian classification, oxisols are. 
generally classified as latossolos amaxelos, alicos, and textura argilosa. In the 
USDA Soil Taxonomy, oxisols correspond to the aplic acrorthox. From soil 
moisture measurements, the porosity was estimated as 0=:O.64, the wilting point. 
as 0 W =O.34, and the field capacity as # r° .51 (Hodnett et.al. 1993). As before, 
the critical point was taken as 70% of the field capacity (0 cr =O.36). Based on 
field work by the Institute of Hydrology (U.K.) and INPE (Rocha 1993, personal, 
communication), the saturated hydraulic conductivity is K s =*79.2 mm/h, the 
pore size distribution index is Bp=0.053, the root depths of upper and lower 
layers are f^[l]=1.0 and f 2 (l]= 0 . 0 , and the soil. thermal conductivity, and soil, 
heat capacity are 0.69 Wm’*K**-and 7.74 x 10^ Jm’^K'V respectively. According 
to Wright et al. (1992), the displacement height for the entire area is 0.17 m 
and the roughness length is 0.026 m, and thus dQ[l]=dQ[2]=0.17 m and 
Zq[1]=zq[ 2]=0.026 m. _ 

Since hourly net radiation was measured at the site, but not long-wave 
radiation, albedo and solar radiation were not used in this model application. 
The study by McWilliam et al. (1993) indicated that the leaf area index (LAI) at_ 
the clearing site is about 1.5. In addition, due to the climatic conditions in 
central Amazonia, the LAI is relatively constant throughout the year. Monteith 
and UnsWorth (1990) suggest that for crops during the growing season, the 
minimum stomatal resistance varies from 50 s/m to 100 s/m. The growing 
season minimum stomatal resistance Was measured to be 50 s/m (Bougeault 
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iSvl). Sellers and Dorman (1987) reported that the minimum stomatal 
resistance is about 50 s/m for short vegetation. Following the approach used by- 
Smith et al. (1993) to determine the minimum stomatal resistance at FIFE , site, 
an optimization on the stomatal resistance within the range (50, 100) was 

carried out using the data of the selected 17 days in Mission 1. The optimization 
indicated that r m »^ [l]=50 s/m gave the best results, although the difference in 
the latent heat flux due to r mifi [l]=50 s/m and r min [l]=100 s/m,, for example, 
is small. The absolute relative changes of the mean latent heat and sensible heat 
fluxes during this period in Mission 1 were 8.0% and 12.1% respectively, while 
the relative change in minimum stomatal resistance was. 100%. Thus, the 
stomatal resistance was taken as 50 s/m. The optimization on mmimiim 
stomatal resistance r mil Jl] was not carried out for the 42 days in Mission 2. As 
for the application at FIFE site, the architectural resistance of grass was taken 
as rQ[l]=2.0 s/m (Ducbudre et al. 1993). 

The soil temperatures were measured at four different depths as described 
in Section 3.3.1, and are shown in Figs. 3.11 (from September 25 to October 5.in 
1990), 3.12 (from October 15 to November 1 in 1990), and 3.13 (from June 30 
to September 10- in 1991). From Fig. 3.11, it can be seen that the soil 
temperature, at 0.05 m depth dropped below 0°C around September 30. This 
Was. due to a measurement problem (Wright 1993, personal communication). 
Fig. 3.12 shows the fluctuation of the soil temperatures at 0.05 m and 0.4 m 
depths is suspiciously large from -October 17 to 28. -During Mission 1, the 
variation of the Soil temperature from its mean at 0.4 m depth is about 14 times 
larger than .the variation of the temperature at depth 0.2 m, and . about .4 times 
larger than the variation at depth 0.1 m. Therefore, the temperature 
measurements at the 0.4 m depth in Mission 1 appear inconsistent, so are the 
measurements at 0.05 m depth. The soil temperature measurements at depths 
0.05 m and 0.4 m in Mission 2 indicated alarming drifts. Thus, the Mission 2 
soil temperature at depth 0.4 m could not be used. Since the fluctuation of the 
soil temperature at 0.2 m depth is small (mean 26.8°C and Standard deviation of 
0.6°C .in Mission 1; 27.1°C mean and 0.6°C standard deviation in Mission 2), 
the mean of the soil temperature at the depth 0.2 m was assumed to be a 
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constant and T2=300°K was used in the model. Thus, depth D 2 that 
corresponds to the constant temperature T 2 is 0.15 m. 

Since there Were no streamflow data available in the studied area, the 
hydrological parameters bj, D s , W s could not be estimated through the 
standard hydrological calibration method, nor could the maximum subsurface 
flow D m be estimated by analyzing the dry period streamflow as was done in the 
FIFE application. Instead, Djn was estimated by multiplying the saturated 
hydraulic conductivity by the average soil slope. At the clearing site, the slope 
is about 3°. As for the three hydrological . parameters b-, D s , and W$, they 
were assigned values based on the ranges suggested by Dumenil and Todini 
(1992). The effect of specifying the three parameters in this. way is insignificant 
to the model predicted latent heat and sensible heat fluxes (see Chapter 4). The 
value of jeach parameter used at the clearing site is given in Table 3.3. 

As noted in Section -3u3.1, soil moisture Was measured at 11 depths, and 
thus the average soil water content could be derived for both the upper layer and 
the lower layer. Since the soil moisture was measured on October 16 of 1990 and 
on July 5 of 1991, the average water contents on each day were derived as 
0=0.368 (0.368 m) and 0=0.436 (3.924 m) for the upper and lower layers on 
October 16, respectively, and as 0=0.420 (0.420 m) and 0=0.466 (4.l9i m) for 
layers 1 and 2, on July 5, respectively. These values were used as the initial 
soil moisture contents of the two test runs for Missions l.and 2 respectively. The 
initial value of the intercepted water by the canopy was estimated based on 
precipitation information . for. the days prior to October 16 and July 5. On 
October 14 and 15, the two days before the beginning of the Mission 1 
calculations, there was no rainfall. Therefore, the initial intercepted water for 
Mission 1 was set to zero. For the Mission 2 application, there was Some rainfall 
on July 3 and 4, but there was no rainfall for the four hours before the beginning 
calculation time at 0000 LMT on July 5. Since the maximum intercepted water 
by the given LAI was 0.17 mm and the average evaporation rate at that time 
was 0.18 mm/hr (Hodnett et al. 1993), the small intermittent rainfall with rates 
less than 0.2 mm/hr before 2000 LMT on July 4 should have evaporated within 
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the four hours before 0000 LMT on July 185, and the initial intercepted Water 
for Mission 2 was also set to zero. 

3.4. Model validation at ABRACOS site 

With the model parameters determined above, the two-layer VIC model 
was applied to the hourly data of the 1? days in Mission 1 and the 42 days in 
Mission 2. The model predictions of the hourly latent heat, sensible heat, and 
ground heat fluxes were tested against the hourly observations. Since surface 
temperature was not measured, the model-predicted effective Surface 
temperature was not tested. 

Fig. 3.15 shows the model predicted (dotted) latent, sensible, and 
ground heat fluxes .and the corresponding observed quantities (solid line) for 
Mission 1. During the 17 days of Mission 1, there was little rainfall, except for 
the one on October 24 which reached about 3.6 mm for one hour (see Fig. 3.14). 
During the first 6 days, the observed latent heat fluxes were about 300 Wm* 2 , 
and the sensible, n eat fluxes were about 200 Win’ 2 , while the latent heat and 
sensible heat fluxes both had smaller values on the last 6 days. The model, 
predicts, latent heat and sensible heat fluxes well for the first 12 days. For the. 
remaining 5 days of Mission 1, the model underpredicts the latent heat fluxes 
and overpredicts the sensible heat fluxes on the same days. 

The reason for the underpredictions of latent , heat fluxes at the end .of 
Mission 1 might be due to underestimation by .the model of the soil moisture 
during that time (Fig. 3.19). In fact, for the 17 days of Mission 1, the model 
indicates that* reduction of upper, zone soil moisture (first 1 m depth) is 
accompanied by a reduction of latent heat and an increase in sensible heat, but 
the observed data do not support this.. From Fig. 3.19, it can be seen that the 
observed upper zone soil moisture increased during October 25 to 27 with no 
rainfall events observed during that period.. This increase of soil moisture makes 
it possible to evaporate more during the next few days as is shown in the 
observed latent heat fluxes.. In contrast, the model shows a decrease in upper 
zone soil moisture that is consistent with the no rainfall period. This decrease of 
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soil, moisture reduces the latent heat flux and increases the sensible heat flux as' 
shown in Fig. 3; 15. The increase of observed soil moisture during the no rainfall 
period might be due to moisture movement from the lower zone. Such capillary 
movement is not represented in the current two-layer VIC model. The largest 
difference between model-simulated and observed soil moisture in Mission 1 was 
about 20 mm. 

For the ground heat fluxes, the model predicts the peaks fairly well for 
the first 7 days, and overpredicts the peaks for the remainder of Mission 1. The 
model-simulated ground heat fluxes tend to give larger negative values on the 
days where small. spikes occurred on either observed latent heat flux or sensible 
heat flux such as on October. 17,. 18, 19, 20, 21, etc., while reasonable 
negative ground heat fluxes were predicted on October 22, 23, and 24,. for . 
example, where such small spikes were not present in the measurements. The — 
means and their standard deviations of the model predicted fluxes and the ones 
of the corresponding observed quantities for Mission 1 jure listed in -Table 3.4a. 
From Table 3.4a, it can be seen that the differences between model predicted 
statistics (mean and standard deviation) and the observed statistics are small.. 

The comparison between model predicted and the observed fluxes for the . 
42 days of Mission 2 is shown in Figs. 3.16 (from July 5 to 18), 3.17 (from July 
19 to August 1), and 3.18 (from August 2 to 15). The soil was moister during 
Mission 2 than it. was during Mission 1. The upper zone soil moisture was above 
the critical point Tor the duration of Mission 2, while it was only slightly above 
the wilting point at the end of Mission 1. The observed latent heat fluxes Were 
about 400 Wm' 2 on the days with no rainfall. From Fig. 3.16, it can be~seen 
that the model predicts the latent heat fluxes very well on both raining days and 
non-raining days. In general, the model overpredicts the sensible heat fluxes 
slightly and underpredicts the ground heat fluxes slightly .as well, but the overall 
prediction is quite good during these days. The model predicts the latent heat, 
sensible heat, and ground heat fluxes quite well from July 19 to August 1 (Fig. 
3.17), and from August 2 to. 15 (Fig. 3.18), The same statistics used to 
summarize the model performance for Mission i were commuted for the 42 days 
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in Mission 2 and are listed in Table 3.4b. The large negative simulated ground 
heat fluxes predicted in Mission 1 did not occur in Mission 2 (see figs. 3.16, 
3.17, and 3.18). A closer look at the observed fluxes during Mission 2 indicates 
that the small spikes in the observed latent heat and sensible heat fluxes were 
not present in Mission 2, suggesting that some of the apparent Mission 1 model 
error may in fact have been attributable to measurement errors. 

During Mission 2, the model simulated moistures were larger than the 
observed ones before August 7, and smaller than the observed ones after August 
7. Nevertheless, the largest difference was less than 30 mm. (Pig. 3.10). The- 
statistics of the. soil moisture predictions and observations for Missions 1 and 2 
are given in Table 3.5. 

3.5. Summary of the model application 

The two-layer VIC model performed well for a hot and. moist condition, 
which is much different from the FIFE site.' The applications to the.ABRACOS 
PIPE sites suggest that the two-layer VIC approach, coupled with a 
simplified vegetation model, may be sufficient to represent land surface fluxes. 
Nonetheless, it must be emphasized that the model testing to date is for two 
small areas with the specific land covers, and climates; -more testing will be 
required at other sites where detailed surface flux data are available before the 
model can be considered to be globally validated. This latter concern, however, 
Ad not limited to this, model alone; a major thrust of such projects as GCIP 
(GEWEX (global energy and water cycle experiment) continental-scale 
international project), and large scale field experiments such as BOREAS 
(boreal ecosystem-atmosphere study), is to provide better large area surface 
moisture and energy flux data for validation of GCM land surface algorithms. 
The approach described in Chapter 2 may be considered as a candidate protocol 
for future validations of GCM land surface parameterizations. 
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Table 3.L Average monthly NDVIs at the FIFE site 


Mo-JF MAMJ J AS OND 


NDVI 53 58 66 89 132 147 145 136 . 122 84 66 62 


Table 3.2. Model parameters at the FIFE site 


Parameter 

Value 

Parameter 

Value 

b i 

0.008 

r Q [l] (s/m) 

2.0 

D m (mxn/h) 

0.34_ 

(*M) 

100.0 

D s .... 

7.7 xlO* 5 

do[l] (m) 

0.25 

W s 

0.96 

2 o[l] (m) _ _ 

0.07 

Bp 

0.16 

C v [l] 

1.0 

K s (mm/h) 

6.44 

fill] 

1.0 

Wf (mm) 

250.0 

m 

0.0 

W£ (mm) 

1250.0 

T 2 ( °K) 

293.6 

Wj 7 (mm) 

G.lSWj* 

a[l] 

0.2 

W[ r (mm) 

0.46W[ 

*[1] (Wm-lK’ 1 ) 

0.514 

$ 

0.5 

C s [l] (Jm-SK- 1 ) 

2.13 xlO 6 
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Table 3.3. Model parameters at the AERACOS site 


Parameter 

Value 

Parameter 

Value 

b ; 

0.1 

r 0 [l] (s/m) 

2.0 

D m (mm/h) 

4.15 

r minM ( s / m ) 

50.0 

Ds 

0.008 

doW (m) 

0.17 

W s 

0.9 

z 0 [l] (m) 

0.026 

B P 

0.053 

Cv[l] 

0.84 

6 

0.64 

LAI[1, m] 

1.5 a 

K s (mm/h) 

79.2 

flW 

1.0 

(mm) 

640.0 

yi] 

0.0 

W '2 (mm) 

5760.0 

t 2 (°K) 

300.0 

Wj v (mm) - 

0.53W? 

*[1] (Wm^K* 1 ) 

0.69 

Wf (mm) 

0.56W? 

C s [l] (Jm^K* 1 ) 

7.44 xlO 5 
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Table 3.4c. Comparison of statistics for. ABRACOS Mission 1 


Names . 

Mean(cal) 

(Wm* 2 ) 

Std(cal) 

(Wm* 2 ) 

Mean(obs) 

(Wm* 2 ) 

Std(obs) 

(Wm' 2 ) 

Latent 

67.63 

98.88 

64.46 

102.16 

Sensible 

40.38 

69.91 

41.19 

67.99 

Ground 

-0.33 

27.S0 

2.07 

18.52 


Table 3.4b. Comparison of statistics for ABRACOS Mission 2 


Names 

Mean(cai) 

(Wm* 2 ) 

Std(cal) 

(Wm* 2 ) 

Mean(obs) 

(Wm* 2 ) 

Std(obs) 

(Wm* 2 ) 

Latent 

93.73 

128.86 * 

95.86 

134.69 

Sensible 

30.98 

44.87 

24.35 

38.31 

Ground 

-4.05 

21.63 

0.43 

22.41 
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Table 15. Comparison of upper zone soil moisture for 
ABRACOS Missions 1 and 2 


Names 

Meanfcal) 

(mm) 

Std(cal) 

(mm) 

Mean(obs) . 
(mm) 

Std(obs) 

(mm) 

Mission 1 

343.34 

9.17 

354.61 

8.78 

Mission 2 

402.64 

• 19.10 

395.06 

14.49 



60 



• MAM 

• 6-ocp 

ti> SP-Suc«t pam 
® SO-Swotr OCP 

• 6-6ow«n 8*tn M* **M r*rh*n > 

X t-EMy CgmWgn hta Mwmwiim 
— Sit* lauridary 
•mill 1-70 
H-177 


Fig. 3.1 Schematic representation of -the FIFE site -with approximate 
locations of King's Creek catchment (shaded area), fluxes stations, 
and meteorological stations (After Famiglletti et al. 1992). 
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Fig. 3.4 Comparison of predicted (dotted) and observed (solid) surface fluxes 
and surface temperature at the FIFE site for June 30-July 11, 1987 
(IFC 2). 
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Fig. 3.5 Comparison of predicted (dotted) and observed (solid) surface fluxes 
and surface temperature at the FIFE site for August 9-20, 1987 
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Fig. 3.C Comparison of predicted (dotted) and observed (solid) surface fluxes 
and surface temperature at the FIFE site for period of extremely 
dry soil moisture during October 6-16, 1987 (IFC 4). 
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Fig. 3.7 Schematic representation of the time frame for the Anglo-Brazilian 
Amazonian Climate Observation Study (ABRACOS) (After 
Shuttleworth et al. 1991). 
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Fig. 3.9 Vertical schematic representation (not to scale) of the location of 
the profile tower (9 m in height) and the eight neutron probe tubes 
(over about 500 m horizontally) used for soil water measurements at - 
the ABRACOS sice (After Hodnett et al. 1993). 
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Fig. 3.13 Soil temperatures at different depths for ABRACOS Mission 
(June 30-September 10, 1091). 
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Fig. 3.14 Hourly precipitation for ABRACOS Mission 1 (October. 16- 
November 1) and Mission 2 (June 30- August 15, 1991). 
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at the ABRACOS site for October 16*NoVember 1, 1990 portion of 


Mission 1. 
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Fig. 3.16 Comparison of predicted (dotted) and observed (solid) surface fluxes 
at *he ABRACOS site for July 5-18, 1991 portion of Mission 2. 
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Fig. 3.17 Comparison of predicted (dotted) and observed (solid) surface fluxes 
at the ABRACOS site for July 19 — August 1, 1991 portion of 
Mission 2. 
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Fig. 3. IS Comparison of predicted (dotted) and observed (solid) surface fluxes 
at the ABRACOS site for Augus t 2-15, 1991 portion of Mission 2. 
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rig. 3.19 Comparison of predicted and observed upper zone soil moisture at 
the ABRACOS site for October 16-November 1,- 1990 portion of 
Mission 1 and for July 5- August 15, 1991 portion of Mission 2. 




CHAPTER 4 SENSITIVITY ANALYSIS 

In Chapter 3, the two-layer version of tLe variable ir. nitration capacity 
(VIC) model (described in Chapter 2) was evaluated by comparing model- 
simulated surface heat fluxes with those observed at two sites: a native tallgrass 
prairie of the FIFE site in the United States ..and a deforested tropical forest of 
the ABRACOS site in Brazil. However, further validation opportunities are 
limited due to the small number of locations globally where surface fluxes of 
moisture and energy have been measured. 

An alternative approach to model testing is sensitivity analysis. Althou gh 
lacking the benefit of comparison with observations, a systematic analysis of the 
mode! sensitivity can at least help us understand how the parameters affect the 
model results. The most commonly used sensitivity analysis method is the one- 
factor-at-a-time approach. The other,, recently suggested for climate model 
assessments by Henderson-Sellers (1992, 1993), is the so called factorial design 
method (Box. et al. 1978) which has the advantage of. investigating *nd 
identifying multiple factor interactions among the parameters. In this chapter, 
the two methods are combined, to investigate and identify the model' parameters 
that most strongly affect the two-layer VIC model. 

4.1. Factorial designs 

Unlike the "change one-factor-at-a-time” approach, the factorial design 
method tests both the sensitivity to changes in individual, parameters, and to 
interactions between groups of parameters. . A general factorial design tests a 
fixed number of possible values for each of the model parameters, and then 
investigates and identifies the ranks of effects of each parameter by running the 
model through all the possible combinations of the parameters. For example, if 
there are k parameters (p^, P 2 , ..., p k ) in the model, and there are 1^ possible 
options (called levels) for the first variable (p x ), 1 2 options .for p 2 , .... and J k 
options for p k , then, a complete factorial design Would include 1^ xl 2 x • • • xl k 
combinations. In this analysis, however, a factorial design at only two levels 
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was considered. Thus, there will he only 2 k combinations of the model 
parameters in the above example. 

Consider a k=3 parameter (2^ factorial) design as an example to illustrate 
the method. Assume that the three parameters are roughness length z Q , leaf 
area index LAI, and the minimum stomata! resistance r^. If the prediction 
variable of interest is annual latent heat flux, the design matrix is: 

z 0 kAL r n - n /> w LgE 

h 


+ ’ + - y4 

~ ~ . + y 5 

+ - + y6 

+ + Y? 

+ + + y g 

where ”+” and ” signs represent the two possible values of each parameter, 
with ”+”-for high values and ” for low values. With this design matrix, the 
effects due to each parameter alone can be estimated as, 




n 


X<Vn> 

i = l 

FT 


(4.i) 


where Ej represents the effect of jth factor (i.e., in jth column), n is the total 
number of experimental runs (i.e., n=3), S~ represents the sign in. row i and 
column j, yj represents the value (e.g., annual latent heat flux) obtained from 
the ith experimental run, and Nj is the number of signs in column j. Using 
Eq. (4.1) and the above design matrix, the effects of parameter interactions on 
the model results can also be estimated based on the signs for the parameter 
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interactions as shown below. The signs are obtained by using the rule that plus 
(”+”) times minus (”-*) gives a minus (”-”), and that minus times 

minus or plus (”+”) times plus (”+”) gives a plus (”+”). 


Run 


1 

2 

3 

4 

5 

6 

7 

8 


Zq • LAI LAI-r^ 

+ + + 

+ 

- + - 

+ - - 

+ 

— — 1 - — 

- - + 

+ + + 


z n -LAI-r . 
u . mm 


+ 

+ 


+ 


+ 


With all the Ejs estimated from Eq. (4.1), the degree of importance of the 
parameters and their interactions can be determined. One way of doing this is to 
plot the Ejs on a normal probability Scale (Box et al. 1978). Any points that are 
outliers from the straight line on normal probability paper could be considered to 
affect the model results significantly, since it is assumed that changes in levels 
of the variables have no real effect on the model results, and thus the model 
results (i.e., Ejs) would occur simply as a result of random variation, about a 
fixed mean. The others are noise, assuming that all the higher interactions are 
negligible. Neglecting higher order parameter interactions is conceptually similar 
to neglecting higher order terms in a Taylor expansion (Box et al. 1978) 
Another way of identifying the parameter's which have major effects on the 
model results was suggested by Henderson- Sellers (1992, 1993). She used an 
iterative method to find thresholds that were two, three, or four standard 
deviations from zero. Any Ejs greater than the estimated thresholds were 
considered to have significant effects on the model results. 


4.2. Fractional factorials 

From the experimental design described above, it can be seen that the 
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number of experimental runs increases geometrically with the number of model 
parameters. A land-surface parameterization scheme with 10 parameters, for 
example, would need 2 10 =-.1024 experimental runs to investigate all the single 
and interacting parameter effects. This number of runs is computer intensive, if 
not totally impossible. However, from a mathematical point of view, it is not 
necessary to conduct all of the experimental runs, since not all of the parameter 
interactions would have appreciable effects on the model results. In fact, there 
tends to be a hierarchy in terms of the magnitude of the parameter effects. The 
single parameter effects (called main effects) tend to Lave greater absolute 
magnitudes than the. two-parameter, interactions, and the two-parameter, 
interactions tend to have greater effects than the three-parameter effects, and so 
on. The terms of a Taylor, series expansion of a response function can be 
associated with the main effects and the interactions (Box.et al. 1978). 
Therefore, the higher orders (or higher interactions) in the Taylor expansion 
senes have larger opportunities to be ignored. In addition, when a large number 
of experimental, runs is introduced, some of the runs are redundant (Box et al. 

19,8' Therefore, only a fraction of the experimental runs of the factorial design 
is needed. 

When a fractional factorial design is used, it implies a tradeoff between, 
the loss of information about higher order interactions and the number of 
experimental runs. In addition, it will introduce confounding patterns where 
certain combinations of the parameters are indistinguishable from others. These 
confounding patterns can be between single parameter and . two parameter 
interactions, two parameter interactions and other two parameter interactions, 
two- and three-parameter interactions, and. so on, . depending on the 
"resolution” at which the fractional factorial experiment is designed. A design 
with resolution R is defined as the one in which no k-parameter effect is 
confounded with any other effect containing fewer than R-k parameters (Box et 
al. 1978). For example, a design with resolution 4 won’t have any main effects 
confounded with any two parameter interactions. However, the main effects can 
be confounded with three parameter interactions, the two parameter 

interactions with other two parameter interactions, and so on. When 
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confounding patterns occur, it is usually considered that lower order effects are 
more likely than higher order effects. However, to be sure that this is the case, 
either further experiments need to be conducted, or physical reasoning is used to 
eliminate certain parameter(s) within the confounding sets. 

The procedure for conducting a fractional factorial design is: 

(1) . Determine the number of model parameters that need to be investigated 
and the number of experimental tuns to be conducted; 

(2) . Determine the design resolution; 

(3) . Select the parameters to which signs will be assigned ("primary 
parameters"), and specify their plus and minus signs; 

(4) . Determine the signs of the remaining parameters that are selected, in (1) 
based on defining relations (defining relations are the equations through, which 
the plus and/or minus signs of the rest of the parameters are determined based 
on the signs assigned to the primary parameters); 

(5) . Write out the design matrix, With ”+” and ” - ” signs representing the two 
possible levels of each parameter; 

(6) . Calculate the effects (also called contrasts) by using Eq. 4.1; 

(7) . Find the confounding patterns, and check if further experiments are 
needed; 

(8) . Rank the degree of importance of these investigated parameters based on 
their absolute magnitudes of effects. 

4.3. Fractional factorial experiments with the two-layer VIC model 
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The two level fractional factorial experiment method is used in this 
Section to investigate the sensitivities of the. parameters in the two-layer VIC 
model. An off-line mode (i.e., prescribe the input forcing data such as solar and 
long-wave radiation, precipitation, Wind, air temperature, and specific 
humidity) sensitivity testing was used. Off-line testing has been reported to be 
an efficient way of exploring the effects of model parameters for GCM land- 
surface parameterizations (Dickinson and Henderson-Sellers 1988). Henderson- 
Sellers (1992, 1993) used the second year results from a two year simulation 
period to analyze the parameter sensitivity of BATS. Use of a Second year 
simulation following a Warm up” year is preferable to the shorter period used in 
some previous studies, e.g., 10 days (Wilsor et al. 1987) or 150 days (Dickinson 
and Henderson-Sellers 1988) since initialization effects are removed (HenderSon- 
Sellers 1992, 1993). 

In this analysis, the one. year of PILPS (Project for Intercomparison of 
Land-surface Parameterization Schemes) prescribed atmospheric forcing 
representing two climatic regimes Was used: . (1) moist, tropical forest; and (2) 
a midlatitude grassland. The PILPS atmospheric forcing data Were obtained 
from the NCAR CCMl-Oz (Henderson-Sellers et al. 1993) for a forested grid cell 
centered at 3°S, 60°W, and a grassland grid cell centered at 52°N, 0°E 

respectively. _ The atmospheric forcing data include (1) downward shortwave 
radiation, (2) downward long-Wave radiation, (3) precipitation, (4) air 
temperature, (5) wind speed, (6) surface pressure, and (7) specific humidity. 
Pitman et al. (1993) provide details of the PILPS experiments , where the above 
two sets of forcing data (i.e., at the tropical forest and midlatitude grassland) 
were used by twenty different land-surface .schemes run to equilibrium. The 
simulated latent and sensible heat fluxes, surface temperature, runoff, and 
other surface fluxes and state variables predicted by the twenty schemes are to 
be compared. The tWo-layer VIC model is among the twenty participants. 
From the monthly precipitation forcing data shown in Fig. 4.1, it can be seen 
that there is considerable seasonality in the precipitation at both sites. The 
monthly precipitation is highest in February, October and December, and 
lowest in May-August at the forest site. For the grassland site, monthly 
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precipitation is highest in July-September, and uniformly lower in the rest 
months. . 

There are 23 parameters associated with the canopy/land-cover, types in 
the two-layer VIC model (Table. 4»1.). In the PILPS study, parameters for 
Manabe bucket, SiB-, and BATS- type models were suggested for both 
grassland and forest sites. Although the two- layer VIC model is not strictly 
similar to any of these types, there are some common, parameters. Therefore, 
the VIC model parameters were determined as was most appropriate based* on 
the given information. For example, the model pore size distribution index Bp, 
and the model maximum subsurface flow D max can be expressed in terms of the 
PILPS parameters as, 


B P = ir 
Dmax = Kg • tan* 

where B is the soil Wetness exponent, and tanx is the surface slope. In addition, 
the soil moisture critical point was taken as 70% of field capacity, and the 
surface albedo (snow-free) was obtained by weighting the albedos with the 
corresponding fractional coverages of vegetation, and bare soil. There, remained 
four parameters that could not be determined based on the PILPS information. 
They are the architectural resistance r 0 , the infiltration shape parameter b-, the 
fraction of the maximum subsurface flow D s , and the fraction of the maximum 
soil moisture content in the lower layer W s . The values of the architectural 
resistance can be obtained from the literature. They are -25 s/m. and 2 s/m for 
the rain forest and grassland respectively (Ducoudre-et al. 1993). The values for 
the other three parameters, however, cannot be estimated well unless 
streamflow data are available and are via a hydrograph fitting procedure. Since 
no streamflow data were available, their values were specified according to 
Dumenil and Todini (1992). The values of the twenty-three VIC model 
parameters jure listed in Table 4.1. 
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Of the 23 parameters, 12 are either relatively easily estimated from field 
data or have relatively narrow feasible ranges compared with the remaining 
parameters. Therefore, these 12 parameters were fixed at the nominal values 
given m Table 4.1. The two level fractional factorial sensitivity analysis was 
only applied to the remaining 11 parameters, where their two levels (high and 
low) were determined based on 50% perturbation about the values given in Table 

4.1. Therefore, the high level ”+” is related to the value given in Table 4.1 for 
the parameter plus 50%, and the low level ” corresponds to the value given 
m Table 4.1 minus 50%. If the values (high and low) determined in this way 
exceeded either the range of that parameter reported in the literature or its 
physical range, then the high and/or low values were adjusted to the literature 
or physical bound. For. the three parameters (b-, D s , and W g ), the two levels 
were obtained based on the ranges suggested by Dumenil and Todini (1992). 
Adding or subtracting 50% error to the parameters may be an exaggeration of 
the parameter estimation errors in some cases, but their ranges should be 
helpful in exploring the parameters sensitivity of the two-layer VIC model. The 

eleven parameters with their high and low values are given in Table 

4.2. 

Four metrics were selected for investigating the sensitivity of the two-, 
layer VIC model. They are the annual total evaporation (mm/yr), annually 
averaged sensible heat flux (Wm* 2 ), annual total runoff (mm/yr), and hourly 
minimum surface temperature (K). To eliminate initialization effects, the 
PILPS procedure of running the simulations for n-year repetitions of the one year 
forcing was followed until convergence was reached. In this case, convergence is 
taken to occur when the monthly maximum, minimum, and mean of the land 
surface fluxes, surface temperature, and runoff were essentially identical for the 
subsequent simulation years. The total evaporation and runoff are the hourly 
accumulated quantities for the equilibrium, year; the sensible heat flux is an 
hourly average of the equilibrium year; and the minimum surface temperature is 
the single lowest hourly temperature simulated for. any hour in the equilibrium 
year. These four measures were selected to represent the radiative and 
hydrological characteristics of the two-layer VIC model. 
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Since there are 11 model parameters in the sensitivity analysis, there 
Would be 2 =2048 experiment runs if a full factorial experiment were used. 
However, this number of funs is fiat necessary as discussed earlier. A fractional 
factorial experiment of 32 runs with resolution 4 (i.e., 2 I y‘ 6 = 32) Was designed 
in this analysis. Using the notation of Box et al. (1978), in the 2j V design, 
there are 5 primary parameters, and the remaining 6 parameters are associated 
with them. The five primary parameters were selected arbitrarily as roughness 
length Zq, soil moisture content at critical point 0 cr , minimum stomatal 
resistance r m j n ) soil thermal conductivity *, and LAI. Let us represent the 11 
parameters by parameter index numbers where the symbol 10, for example, 
represents the tenth parameter, then we can specify them as (see Table 4.2), 

z 0 " 1 *cr - 2 r^ - 3 * - 4 LAI - 5 

c s- 6 0-7 0* - 8 bj - 9 D s - TO W s - IT 

The last 6 parameters are related to the first 5 primary parameters through the 
defining relations given by Box et al. (1978). The defining relations are: 

I = 1236, I = 2347, I = 3458, (4.2a) 

I = 1349, I = 14510* I = 245H. . (4.2b) 

where I represents a column of all plus signs, and 1236 represents the resultant 
sign from the signs of parameters 1, 2, 3, and 6 based on the plus times minus 
rule discussed earlier. In other words, the above equalities mean that the signs 
on both sides are the same. For example, if the Signs for 1, 2, 3, and 6 are 
and ” respectively for run 1, the sign of 1236 is then ”+” 
according to the plus times minus rule, and it . has the same sign as T.. If the 
sign of 1236 is then it won’t be equal to T, i.e., 1^1236. From_the 

above six defining relations, the sign of a parameter index number can be 
determined if the signs of the remaining variables are known. This can be done 
by multiplying the number (i.e., parameter) on both sides of an equation. For 
instance, the sign of number 6 can be determined by 6=1 -6=1236 -6=123 1=123. 
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Since the first five numbers have been chosen as the primary parameters, their 
signs can be specified, and the signs for the remaining numbers can be then 
obtained by multiplying its number on both sides of the equation. The six 
relations from the above defining relations are, 

6 = 123, 7 = 234, 8 = 345, (4.3a) 

9 = 134, IU = 145, U = 245. (4.3b) 

The above six relations are called generators. Thus, by specifying the plus and 
minus signs to the five primary parameters for the 32 experimental runs, the 
signs of the remaining six parameters in the 32 runs can be obtained through the 
above six generators. The relations in the above six generators are not unique. 
There are other ways to define the six generators for the il parameters. 
However, the design shown here is preferred because it ensures that all the main 
effects are not confounded with the effects of any of the two parameter 
interactions. Only two parameter interactions are confounded with each other if 
higher order effects are not considered (Box et al. 1978). In other words, the 
design shown here guarantees a design with resolution IV in 32 experimental 
runs. The design matrix wits all the plus and minus ” — ” signs for__each of 
the 11 parameters in the 32 runs, is shown in Table 4.3. 

From the above six defining relations, all the confounding patterns of this 
resolution IV design with 32 experimental runs for 11 parameters can be found 
based on the procedures described by Box et al. (1978). In the following, 
however, only the confounding patterns for the two parameter interactions are 
discussed; higher order interactions . are assumed to be negligible. Multiplying 
two of the six defining relations, at a time (e.g., I = I- 1 = 1236*2347 = 

1 *2 *2*3 *3 <467 — 1*1*1*467 = 1467) gives the I’s With only four parameter 
combinations as follows, 

I = 1467 = 2469 = 2578 = 1279 = 357U 
= 1589 = 13810 = 238U = 359l0 * 12lU Tl (4.4) 
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Multiplying three at a time gives the I’s with only four parameter combinations 
as follows, 


I = 3679 = 56710 = 47811 = 489lU = 3610 IT 
= 26810 = 168H = 569IT (4.5) 

Multiplying four at a times gives the I’s with only four parameter combinations 
as follows, 

I = 791011 (4. 6 ) 

Multiplying five and six at a time won’t end up with any I’s that can be 
expressed by only four parameter combinations. Therefore, they cannot result 
in any. confounding patterns among two parameter interactions. By multiplying 
12 through Eqs. 4.2, 4.4 — 4.6 and omitting words with three or more numbers, 
for example, we obtain 

I * 12 = 1236-12 = 1279-12 = 12I0U-12 

i.e., 

12 = 36 = 79 = 10 lT 


Thus, the two parameter interaction 12 is confounded with interactions of 36, 
79, and 10 11. Following this procedure, the confounding patterns for the two 
parameter, interactions for the entire design can be obtained based on the 
relationships defined in Eqs. 4.2, 4.4 — 4.6. The confo unding patterns for two 
parameter interactions are: 


12 = 36 = 79 = T6 IT 
14 = 39 = 5l0 = 67 
16 = 23 = 47 = 8lT 
18 = 3l0 = 59 = 6IT 
110 = 2U = 38 = 45 
24 = 37 = 511 = 69 


13 = 26 = 49 = 8l0 
15 = 410 = 89 
17 = 29 = 46 
19 = 27 = 34 = 58 
111 = 2l0 = 68 
25 = 4ll = 78 
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28 = 311 = 57 = 610 . 35 = 48 = 711 = 9l0 

56 = 7l0 = 9ll ... ( 4 ^ 

The matrix having all the plus and minus signs for the two parameter 
interactions are shown in Table 4.4. 

4.4. Results of fractional factorial experiments 

With the 2jy design discussed above, the thirty- two experimental runs 
were conducted for the grassland site and the forest site respectively. The results * 
for the four metrics of the 32 two-layer VIC model runs are listed in Tables 4.5a 
and 4.5b for grassland and forest respectively. The eleven main effects and the 
two parameter interaction effects (specified in Eq. 4.7) are listed in Tables 4.6a 
and 4.6b for grassland and forest respectively. The iterative method (Henderson- 
Sellers 1992, 1993) discussed in Section 4.1 Was used to find the two-, and/or . 
three-standard deviations (2a and 3a). All the main effects and the two 

parameter interaction effects that are greater than 3a in any of the four measures 

at both sites are shown in Tables 4.7a and 4.8a respectively. 

From Table 4,7a, it can be seen that only LAI was significant for three 
metrics (the total evaporation and runoff, and the hourly average sensible heat 
flux). For the minimum hourly surface temperature metric, two single 
parameter effects (k and C s ) and a two-parameter interaction effect (/cC s ) were 
detected to be important. In the forest Case (Table 4.8a), on the basis of the 3a 
criterion all the effects were insignificant on all metrics except for the metric of 
minimum hourly , surface temperature where *, C<j, and were found to be 
significant. The effect of xCg are confounded by the effects of zq0 and 0 cr b- as 
shown in Eq. 4.7. In general, additional small full factorial experiment runs can 
be -conducted to resolve these aliases if the physical reasoning is not obvious. In 
this case, however, since. Zq and 0 , and 0 cr and b* are much less correlated than 
k and C s , and also Zq, 0, 0 cr , and bj have insignificant main effects compared 
with that of k and C s , it is likely that the greatest effect in the confounding 
patterns is from «C S rather than from 5jj0 or 0^bj. Thus, the interaction 
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between k and C s was considered to have the greatest effect on the minimum 
surface temperature metric. 

At the 2 a threshold, more parameters were shown to have large effects on 
the four Selected metrics (Tables 4.7b and 4.8b). From Tables 4.7b and 4.8b, 
the importance of the single parameters and the two-parameter interactions can 
be ranked based on absolute values of their effects. These identified parameters 
were classified into two categories, primary importance and secondary 
importance as shown in Tables 4.9 and 4.10. The effects of the parameters were 
also identified by plotting the values in Tables 4.6a and 4.6b on a normal 
probability scale. The outlier's shown on the plots (Figs. 4.2 and 4.3) are the 
parameters that have significant effects. From Fig. 4.2, it can be seen that for 
all the metrics except for the minimum hourly ..temperature metric, the outliers . 
that are identified by the. probability scale approach are the primary ones 
detected by the threshold, method for the grassland site. For the forest site, all 
the primary and secondary outliers from the threshold method are identified by. 
the probability Scale approach. 

Although the model parameters and their ranks that were identified .are 
not exactly the Same for the grassland and forest Sites (Tables 4.9 and 4.10), 
there is some consistency in the results.. The leaf area index LAI, porosity 0, 
the minimum stomatal resistance r^^, critical point and wilting point of soil . 
moisture have important effects on the first three metrics of the two-layer VIC 
model at both sites. For the grassland site, the roughness length Zq and the 
interaction between porosity 6 and the minimum stomatal resistance r m > Q are of 
secondary importance for the same three metrics. For the minimum hourly 
surface temperature, the soil thermal conductivity and soil heat capacity were of. 
primary- importance at both sites, and the interaction between the two is 
secondarily important. 

Henderson-Sellers (1992, 1993) used, the fractional factorial method to 

analyze the parameter sensitivities of BATS. BATS has 23 model parameters 
that are related to the canopy/lahd-cover types. In the interquartile-range 
ecotype parameter experiment, ten parameters were kept constant, and thus, 
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the experiments had 13 parameters with two different levels each. The metrics 
used Were the maximum canopy temperature, . evapotranspiration, annual 
runoff, the minimum daily upper-layer soil temperature, and the net carbon 
gain. Based on the thresholds 4a and 3 a, the Vegetation roughness length, soil 
porosity, and stomatal resistance to vegetation light sensitivity were found, to be 
primarily important, while vegetation albedo, soil color, wilting point, 
minimum leaf area index, soil moisture diffusivity, and the soil thermal 
conductivity Were found to be of secondary importance (Henderson-Sellers 1992, 
1993). In the full-range ecotype parameter experiment runs with the same 
threshold criterion (Henderson-Sellers 1992), the maximum leaf area index, 
roughness length, shortwave vegetation albedo, the interaction of roughness 
length and albedo, the stem-area index, and the near-infrared vegetation albedo 
were found to be important. Comparing the results described in this analysis 
with the ones obtained by Henderson-Sellers (1992, 1993), it can be seen that 
among the parameters common to both, models, the. leaf area index, porosity,, 
roughness length, soil thermal conductivity, and the parameter related to the 
stomatal resistance are important, despite the significant differences in the 
structure of BATS and two-layer VIC. 

From the study results, it can be seen that the experiments are highly 
fractional. That is, there are many higher order interactions that were not - 
studied. There may exist three or higher order parameter, interactions that 
require examination. Even though this analysis is preliminary, it reveals some 
common parameter effects that are important to the two-layer VIC model, and 
perhaps, to other land-surface schemes as well. 

4.5. One-at-a-time analysis 

From the fractional factorial analysis, it is seen that the infiltration shape 
parameter bj, the fraction parameter of maximum subsurface flow D s , the 
fraction parameter of the lower layer maximum soil moisture. W s , and their 
interactions do not significantly affect the four model metrics (Tables 4.6a and 
4.6b). Among the parameters that show important effects, the above analysis 
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indicates that the ranks of some parameters at the grassland site are different 
from those at the torest site. For example, at the forest site, the critical soil 
moisture had the least influence on the annual total evaporation, and the 
minimum stomatal resistance parameter had the greatest influence, of all the 
parameters tested. At the grassland site, the critical soil moisture affected 
evaporation more than the minimum Stomatal resistance. To explore the 
characteristics of these parameters further, the one-at-a-time method was used 
to allow investigation of more than two levels. 

First, the three parameters bj, D s , and W s were investigated at four * 
levels each. The four different values for each of the three parameters (Table 
4.11) were specified based on the ranges suggested by Dumenil and Todini 
(1992). The annually averaged latent heat flux (Wm'^), sensible heat, flux 
(Wm" 2 ), annually averaged surface temperature (K), annually total evaporation . 
(mm/yr), and the total runoff (mm/yr) from the equilibrium year axe shown in 
Tables 4.11 and 4.12 for both grassland and forest sites respectively. Fig. 4.4 
shows the monthly means of the control and the sensitivity tuns for latent and 
sensible heat fluxes, the surface temperature, and the monthly total runoff at 
both sites. The control runs shown in Fig. 4.4 .Were obtained by running the . 
model to the equilibrium year with the parameter values given jn Table 4.1. 

From Fig. 4.4, and Tables 4.11 and 4.12, it is clear that the infiltration 
shape parameter bj does not Significantly affect latent and sensible heat fluxes, 
and the surface temperature, especially at the grassland site. At the forest site, 
changing b- did not change the monthly latent heat and Sensible heat fluxes 
much, except in August, where increased b- values reduced the amount of 
water infiltrated into the soil, and thus decreased the latent heat flux and 
increased the sensible heat flux. The apparent change in latent heat flux in 
August was because it Was very dry in that month (see Fig. 4.1), so any change 
in the infiltration amount would strongly affect the water available for 
evaporation. In comparison, bj has stronger effect on runoff. Figures 4.$ and 
4.6 partition the runoff into surface and subsurface flow for the range of bj values 
for the grassland and forest sites respectively. The total monthly runoff (Figs. 
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4.5. and 4.6), especially the totally annual runoff (Tables 4.il and 4.12), Were 
much less sensitive to bj. This Was due to the partitioning into surface and 
subsurface flow. D s and W s had minimal effect on the monthly and annual 
latent and sensible heat fluxes, and the surface temperature (Tables 4.11 and 
4.12). Only the monthly runoff distributions vary with different parameter 
values - (Fig. 4.7). Tables 4.11 and 4.12 show that, although the total annual 
runoff is insensitive to Dg and Wg, the monthly runoff distribution varies 
significantly (Fig. 4.7). 

The results from the one-at-a-tim i method confirmed the findings from 
the fractional factorial experiment, and also showed the effect of the three 
hydrologic parameters bj, Dg, and W s on the redistribution of runoff (both 
surface and subsurface). The wetness indices (the .ratios of annual latent heat 
plus sensible heat to annual precipitation), equal to 0.44 at grassland and 0.5 at 
forest, indicated that both sites studied here were climatologically moist. The 
results obtained here might be different from the ones under a dry climate 
condition, especially for parameter bj. 

The one-at-a-time. method Was also used to study the sensitivities of the 
critical point and the minimum stomatal resistance, both of which were found to 
be important in the fractional factorial studies. The results are given in Tables 
4.13 and 4.14, for grassland. and the forest, respectively. Table 4.13 indicates 
that for the grassland site, a change of 7% to 14% in the critical point results in 
almost the same sensitivities as a change of 15% in the mmiirmfn stomatal - 
resistance for all the measures. Table 4.14, however, indicates that at the 
forest site, a smaller change in. the minimum stomatal resistance than in the 
critical point caused larger changes in the latent and sensible heat fluxej, and 
total runoff. Therefore, the sensitivities of the same parameters under different 
climatic conditions can be different. The relative order of importance of the two 
parameters at different sites Was the same as indicated by the fractional factorial 
method. 

4.6. Supplementary fractional factorial experiments 
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From Sections 4.4 and 4.5, it can bo seen that the annual total 
evaporation (mm/yr), annually averaged sensible heat flux annual 

total runoff (mm/yr), and the hourly minimum surface temperature (K) are not 
sensitive to b-, D s , and Wg, but the monthly runoff distributions are Sensitive 
to these parameters. To examine the sensitivities of the three hydrologically 
related parameters further, three different metrics were used. They are the sum 
of the absolute difference of monthly evaporation (mm/mo), runoff (both surface 
and subsurface) (mm/mo), and the sensible heat flux (Wm*^) between the 
control experiment and the 32 experiments described in Section 4.4. 

The results of the 32 experiment rims with the three new metrics are 
shown in Tables 4>15a and 4.15b, for the grassland and forest sites, respectively. 
The parameter effects with the new metrics are shown in. Tables 4.16a and 4.16b, 
for the grassland and the forest cases, respectively. By combining the 3<r and 2& 
threshold approach and the probability Scale approach (Section 4.4), the 
parameters that are sensitive are identified. .From Figs. 4.8 and 4.9 (or Tables 
4.17 and 4.18), it is clear that the effects of parameter interactions axe more 
important than the effect from the single parameters under the new metrics for 
evaporation and sensible heat flux for both the grassland, and forest sites. Since 
the two parameter interactions are confounded with each other in this resolution 
IV experimental design, further, experiment runs would be needed to identify 
the confounding patterns. Rather than performing such runs, the confounding 
patterns that cannot be excluded based on physical grounds are retained as 
possibly important two parameter interactions (see Tables 4.17 and 4.18). 

For the sum of the absolute difference of .monthly runoff (combined 
surface and subsurface) between the control experiment, and the 32 experiments, 
the fractional factorial experiment analysis identified the two parameter 
interactions of bjD s (confounded by 6\V S , r ^ LAl), 0 w bj, and bjWg 
(confounded by 0D S ), and the parameter W s as sensitive at the grassland site. 
For the forest site, the single parameters b- , D s , and W s are sensitive. Like 
the grassland case, the two parameter interactions b-D s (confounded by 0W S , 
r mi nkAI), and b-W s (confounded by 0D S ) are important. These findings are 
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consistent With those obtained from the one-at-a-time analysis discussed in 
Section 4.5. The analysis conducted here indicates that different Metrics used in 
the fractional factorial analysis could result in the identification of different 
parameters and parameter interactions. 

4.7. Additional analysis of model sensitivity experiments 

In Sections 4.3 -4.6, the sensitivity of the two-layer VIC model to the 
model parameters were studied oh the basis of several metrics. In this section, 
the model structure is explored for a few special cases, with the_ model 
parameters fixed at the values given in Table 4.1. The four cases are: 

(1) . Vegetation (grass in grassland, and tree in the forest) covers 100% of the 

•and surface, instead of 80% in grassland and 90% in the forest as used 
previously; 

(2) . As case (1), but without .atmospheric stability correction, and with 
architectural resistance set to zero (see Chapter. 2); 

(3) . As case (1), but the surface was kept wet so that the actual evaporation 
equals the potential evaporation; 

(4) - As case (1) and (2), but with Wet surface. 

The results for the above four caSes-and for the standard case are.given in 
Tables 4.19 and 4.20 for the grassland and forest sites respectively. At the 
grassland site in case (1), the change of the 20% vegetation cover from the bare 
soil to grass reduced the evaporation from 617.21 mm/yr Jp 565.51 mm/yr, and 
increased the annually averaged sensible heat flux -from -5.10 Wm* 2 to 0.02 
Wm , and the total annual runoff from 648.83 mm/yr to 700.54 mm/yr. 
Changes m the annually averaged surface temperature were negligible. At the 
forest site, the change of the 10% vegetation cover from bare soil to trees 
resulted in similar model effects. The annual evaporation decreased from. 1357.59 
mm/yr to 1176.49 mm/yr, the sensible heat increased from 25.93 Wm* 2 to 41.84 
Wm , and the runoff changed from 1909.50 mm/yr to 209Q.55_mm/yr. Again, 
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the. surface temperature changes Were negligible. Elimination of the stability 
correction resulted in only a small increase in evaporation at the grassland site, 
but the effect at the forest site was large. The relative increase of evaporation at 
the forest site was 49.8%, while it Was only 1% at the grassland. Such a big 
difference was partially due to setting the architectural resistance to zero under 
the no stability correction case at both sites. The architectural resistance (see 
Table 4.1) in the forest site (25.0 s/m) was much greater than that at the 
grassland site (2.0 s/m), thus by setting it to zero, it had a greater effect on 
evaporation at the forested site than at the grassland site. In spite of the effects 
of different values of architectural resistance, the influence of the stability 
correction on the latent heat flux Was .apparent, as Was the different effects 
according to climatic conditions. Although the annual total evaporation 
increased at both sites under case (2), the surface temperature at both sites 
increased slightly instead of decreasing. 

By keeping the surface wet in cases (3) and (4), only five of twenty-three 
(see Table 4.1) parameters remained that Would affect the model results. These 
are aerodynamic roughness length, displacement height, surface albedo (snow 
free), soil thermal conductivity, and soil heed capacity. 

By comparing the results (Tables 4.19 and 4.20) of cases (3) and (4), it 
was found again that the effect of the atmospheric stability correction varies 
significantly . with the climatic conditions. For example, the annual total- 
evaporation obtained in case (4) at the forest site was more than three times as 
large as that obtained from case (3). However at the grassland site, the annual 
total evaporation in case (4) was only 1.2 times as large as that in case (3). This 
result implies that under some cases, the atmosphere is closer to the neutral 
condition than the others, and thus, the stability correction is not as important 
for these cases as for the others. 


4.8. Summary of sensitivity analysis 

The studies conducted here show that the combined use of the fractional 
factorial method and the one-at-a-time method is an efficient way of examining 
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the relative importance of model parameters. From the analysis discussed in this 
chapter, it was found that the leaf area index LAI, porosity 6 , and the 
minimum stomatal resistance r min were the most important parameters of the 
two-layer VIC model in terms of the three metrics associated with annual total 
or average surface fluxes, while the roughness length Zq, critical point and 
wilting point of soil moisture, and the interaction between porosity 6 and the 
minimum stomatal resistance r ^^ were of secondary importance for the same 
metrics. For the minimum hourly surface temperature, the soil thermal 
conductivity and soil heat capacity were primarily important, and the 
interaction between the two was secondarily important. 

The studies indicated that the surface fluxes and surface temperature 
Were not sensitive to the fraction of the maximum subsurface flow ana the 
fraction of lower, layer maximum Soil moisture for the four metrics. Both 
parameters resulted in redistribution of the monthly runoff within a year, but 
they did not change the total amount of the annual runoff, tinder relatively wet 
climatic conditions, it was found that the surface fluxes and the surface 
temperature are not sensitive to the infiltration shape, parameter. This 
parameter partitioned the. streamflow into surface . and subsurface flows 
differently based on its different Values, but it did not affect significantly the 
total amount-of amw*l runoff. 

When the sum of. the absolute difference of monthly runoff (both surface 
and subsurface runoff) (mm/mo) from the control experiment was used as the 
metric, the parameters b-, D s , and W s , and their interactions among 
themselves and With parameters 6 or are found to be important and sensitive. 
In addition, when the metrics were the sum of the absolute difference of 
monthly evaporation and sensible heat flux between the control experiment and 
the 32 experiments, similar parameters as those indicated by the metrics of the 
annual total evaporation and annually averaged sensible heat flux were 
identified. However, it was the two parameter interactions of those parameters 
rather than the single parameter effect that were found to be more important 
this time. 
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Finally, the analysis conducted in Section 4.7 showed that the model 
results were sensitive to the inclusion of the stability correction. Thi^ was 
especially the case for the forested site where evaporation was primarily 
atmospherically controlled, and less so for the grassland site where evaporation 
was. mostly soil controlled. 
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Table 4.1. T " 0 -layer VIC model parameters for PUPS grassland and forest sites 


Parameter names 

■ Grass 

Forest' 

Fractional coverage of vegetation Cv 

0.3 

0.9 

Depth of upper layer dj (m) 

1.0 

1.0 

Depth of lower layer (m) 

9.0 

9.0 

Roots in upper layer 

100% 

90% 

Roots in lower layer ^ 

0% 

10 % 

Saturated hydrauhc.conductivity K s (m m/s) 

0.45 xlO* 2 

0.16 xio* 2 

Soil wetness exponent (B parameter) 

6.8 . 

9.2 

Slope x 

0.17 

0.17 

Displacement height dg (m) 

0.0 

18.0 

Surface albedo a (snow free) 

0.21 

0.131 _ 

Constant soil temperature at lm depth T2 (K) 

274.6 

299.6 

Architectural resistance rQ (s/m) 

2.0 

25.0 

Soil porosity d 

0.51 

.0.6 

Fract. of water cent, at which perm. wilt, occurs 

0.37S 

0.487 

Fract. of water cont. at which criti. point occurs 

0.7 

0.6 

Infiltration shape parameter b- 

0.1 

0.03 

Fraction of maximum subsurface flow D s 

0.008 

0.008 

Fraction of lower layer maximum soil moisture W 3 

0.9 

0.9 

Aerodynamic roughness length zq (m) 

0.1 

2.0 

Leaf area index LAI 

Monthly LAIs 

5.0 

Minimum stomatal resistance r ;x ^ a (s/m) 

200.0 

150.0 

Soil thermal conductivity k (Wm'^K**) 

1.03 

0.866 

Soil heat caoacitv C« f Jm* 2 K**> 

2.085x10® 

1.756x10^ 
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Table 4.2a. Range of the eleven parameters at PILPS grassland site 


No. 

Parameter names 

mam 

■nsi 

1 

Aerodynamic roughness length zg (m) 

0.15 

0.05 

2 

Critical point of soil moisture 6 ct 

0.536 

0*179 

3 

Minimum stomatal resistance r„ ; , (s/m) 

200.0 

100.0 

4 

Soil thermal conductivity * (Wm^K^) 

1.545 

0.515 

5 

Leaf area index LAI 

1.5 -LAI . 

0.5 LAI 

6 - 

Soil, heat capacity C s (Jm'%* 1 ) 

3*128 x 10*® 

1.043x10* 

7 

Soil porosity 9 

0.66 

0.33 

s 

Permanent wilting point of soil moisture 0 W 

0.29 

0*097 

9 

Infiltration shape parameter b- 

0.50 

0.10 

10 

Fraction of maximum subsurface flow D a 

0.50 

0.004 

11 

Fraction of lower layer max. soil moist. W a 

0.90 

0.10 
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Table 4 ; 2b. Range of the eleven parameters at PILPS forest site 


No. 

Parameter names 

ffigh(+) 

Lowf -) 

1' 

Aerodynamic roughness length Zq (m) 

3.00 

1.00 

2 

- Critical point of soil moisture 9 a 

0.54 

0.18 

3 

Minimum stomatal resistance r^ ; „ (s/m) 

200.0 

75.0 

4 

Soil thermal conductivity * (Wm^K^) 

1.299 

0.433 

3 

Leaf area' indent LAI 

7.5 

5.0 

6 

Soil heat capacity C s (Jm'^K'*) 

2.634x10^ 

0.878 xlQ* ( 

m * 

Soil porosity 9 

0.66 

0.33 

8 

Permanent wilting point of soil moisture 0 W 

Q.438 

0.146. 

9 

Infiltration shape parameter b- 

0.50 

0.001 

10 

Fraction of maximum subsurface flow D s 

0.50 

0.004 

11 

Fraction of lower layer max. soil moist. W s 

0.90 

0.10 
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Table 4;3. Design matrix of the eleven parameters 


Runs 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

l 
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Table 4.4. 


Computation matrix for two parameter interactions 
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Table 4.5a. Results of 32 rtrns at PUPS grassland site 


Runs 

ET (nvtn/yr) 

SH (Wm-2 ) 

R.(itun/yr) 

Ttnin (K) 

1 

604.56 

-1.56 

661.47 

254.03 

2 

622.51 

-2.59 

643.52 

256.32 

3 

598.10 

-0.67 

667.92 

256.13 

4 

620.12 

-2.80 

645.91 

253.97 

5 

544.37 

2.85 

721.66 

256.15 

6 

566.19 

0.82 

699.24 

254.01 

7 

396.12 

12.10 

869.91 

254.37 

8 

423.29 

10.95 . 

842.74 

256.37 

9 

649.82 

-10.38 . 

616.21 

256.32 

10 

677.21 

-11.11 

588.82 

257.22 

11 

408.62 

6.12 . 

851.41 

257.36 

12 

435.49 

2.89 

830.54 

256.41 

13 - 

531.16. 

• -1.63 

734.87 

251.33 

14 . 

555.19 

-4.62 

710.84 

256.45 

15 

529.11 

-2.84 

136.92 

256.35 

16 

541.37 

-2.31 

724.85 

257.24 

17 

664.66 

-5.54 

601.37 

253.88 

18 

673.62 

-6.00 

592.41 

256.24 

19 

700.76 

-7.33 

565.27 

256.09 

20 

712.10 

-9.16 

553.93 

253.90 

21 ' 

723.20 

-8.74 . 

542.83 

256.10 

22 

762.23 

-12.49 

503.80 

253.93 

23 

585.47 

-0.41 

680 < 56 

254.03 

24 

620.15 

-2.49 

645.88 

256.32 

25 

852.36 

-22.79 

413.87 

256.26 

26 

911.43 

-26.03 

354.60 

257.14 

27 

659.25 

-9.39 

606.18 

257.26 

28 

_ 698.35 

-14.21 

567.18 

256.34 

29 

601.86 

-5.72 

664.16 

257.27 

30 

616.57 

-8.79 

649.46 

•256.39 

31 

626.29 

-9.01 

639.74 

256.33 

32 . 

635.52 

.-8.55 

630.50 

257.22 
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Table 4.5b. Results of 32 runs at PILPS forest site 


Runs 

ET(ntth/yr) 

SH (Wm-2 ) 

R(mm/yr) 

Tmin (K) 

1 

1358.52 

27.73 

1908.55 

278.82 

2 

1333.75 

28.96 

1933.32 

281.91 

3 

1506.33 

17.38 

1760.74 

281.72 

4 

1473.13 

18.78 

1793.94 

278.85 

5 _ 

1149.47 

42.40 

2117.40 

281.90 

6 

1227.79 

36 4 38 

2039.28 

278.91 - 

7 

953.33 

56.26 

2313.74 

279.10 

8 

997.80 

53.30 

2269.27 

-282.00 

9 

1351.37 

24.90 

1915.51. 

284.79 

10 . 

1437.59 

19.64 

1829.48 

285.61 

11 

• 1057.90 

46.42 

2209.17 

286.17 

12 

1201.96 

36.08 

2065.10 

285.60 

13 

1218.36 

35.06 

2048.70 

285.88 

14 . 

1190.92 

36.95. 

2076.14 

285.52 

• 15 

1261.56 

31.13 • 

2005.51 

285.22 

16 

1248.03 

33.46 - 

2020.52 

286.00 

17 

1381.12 

26.10 

1885.94 

278.76 

18 

.1378.51 

25.86 

1888.56 

281.84 

19 

1401.40 

24.74 

1865.67 

281.77 

20 . 

1383.33 

25.12 . 

1883.54 . 

278.93 

21 

1399.09 

25.03 

1867.97 

281.80 

22 . 

1475.02 

18.64 

1792.05 

278.86 

23 

1123.91 

-4-4.18 

2143.16 

279.02 

24 

1232.39 

36.40 

2034.63 

281.93 

25 

1658.31 

3.27 

1610.24 

284.56 

26 

1741.38 . 

-2.41 

1525.69 

285.51 

27 

1363.92 

25.03 

1903.15 

285.82 

28 

1430.62 

19.45 

1836.44 

284.87 

29 

1276.05 

31.08 

1991.02 

285.68 

30 

1250.32 

32.69 

2016.74 

284.84 

31 

1310.30 

27.63 

1956.77 

285.14 

32 

1215.74 

35.76 

2051.13 

286.08 
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Table 4.5a. Parameter effects at PILPS. grassland site 


Index No * 

El> 

(ntin/yr) 

R. 

(min/yr) 

SH 

(Witi-2) 

Tmin 

(K) 

i 

24.30 

-24.80 

-1.97 

0.017 

2 

-35.43 

. 85.43 

5,45 

0.044 

3 

-76.92 

76.92 

4 ,98 

0.058 

4 

7.00 

-6.93 

-5.96 

1.694 

c 

mi 

146.23 

-146.28 

-9.49 

-0.087 

6 

' -0.21 

0.20 

1.01 

1.546 

7 

97.09 

-97.06 

. -6.23 

-0.129 

3 

-67 . 62 

67.60 

4.36 

0.024 

q 

-9.. 73 

9.76 

0.66 

0.023 

l6 

1.73 

-1.75 

-0.15 

-Q.022 

11 

2.57 

-2.59 

-0.17 

0.008 

12 

-1.90 

1.92 

0.19 

-0.028 

13 

-1.35 . 

1.38 

0.21 . 

-0.017 

14 

. 1.35 

-1.35 

-0.16 

-0.018. 

15 

2.23 

-2.31 

-0.38 

0.016 

16 

17.43 

-17.40 

-0.98 

0.031 

17 

0.55 

. -0.55 

0.38 

-0.654 

IS 

-6.23 

6.23 

0 . 47 

-0.011 

19 

-5.07 

5.07 

0.20 

-0.032 

110 

12.99 

-12.99 

-0.58 

0.021 

111 

-0.12 . 

0.12 

0.03 

-0.0i4 

24 .... 

-22 .20. 

22.20 

1.27 

-0.021 

25 

14.49 

-14.51 

-1.01 

-0.009 

23 

3.34 

-3.34 

-0.22. 

0.034 

35 

-10.79 

10.77 

0.55 

0.002 

56 

i.ll .• 

-1.14 

0.01. 

0.027 
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Table 4.6b. .Parameter effects at PILPS forest site 


Far index No. ET 

(riim/yr) 

R 

. (mm/yr) 

etf 

(Wm-2 ) 

Tmin 

(K) 

I 

27.96 

-27.96 

-2.08 

0.070 

2 

-104.12 

104.12 

7.43 

0.190 

3 

-183.07 

183.06 

13.08 . 

0.147 

4 

27.46 

-27.23 

-4.45 

5.074 

5 

123.35 

-123.35 

-9.14 

-0.161 

6 

-4.61 

4.61 

• o.Io 

• 1.865 

/ 

155.65- 

-155.51 

-11.23 

-0.131. 

r> 

3 

-127.58 

127,34 

9.08 

0.052 

Q 

-57.95 

58.08 

4.17 

0.054 

it - 

-16.17 ' 

15.99 

1.10 

-0.040 

11 

11.07 

-11.25 

-0.74 

-0.008 

12 

-2.42 

2.55 

0.28 

-0.031 

13 

-9.72 

9.90 

0.93 

-0.019 

14 

-0.61 

0.61 

0.47 

0.027 

15 

-3.31 

3.57 

0.14 

-0.030 

16 

“1.27 

1.56 

0.06 

-0.051 

17 

-' .44 

7.44 

0.70 

-1.087 

18 - 

-1.45 

1.44 

0.24 • 

-0.006 

19 - 

24.09 

-24.10 

-1.66 

0.033 

lfO 

31.52 

-31.52 

-2.25 

-0.124 

111 

-13.52 

13.52 

0.95 

0.009 

24 

-25,16 

25.16 

1.80 

0.125 

25. 

-33.15 . 

32,92 

2.33 . 

0.025 

23 

-12.19 

12.19 

0.86 

0.119 

35 

1.10 

-1.23 

-0.05 

0.015 

56 - 

4.06 

-4.24 

-0.25 

0.068 
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Table 4.7a. Parameters selected at the PILPS grassland site 

based on a threshold of |3<r| 


Outliers 

Evaporation (mm/yr) 

LAI 

Runoff (mm/yr). 

LAI 

Sensible heat (Wm* 2 ) 

LAI 

Min. surface temperature (K) 

Cs, «C S 


Table 4.7b. Parameters selected at the PILPS grassland site 

based on a threshold of |2<r| 


Outliers 

Evaporation (mm/yr) LAI, 0, 0 CT , r^, 0 W , z Q , r““ 

Runoff (mm/yr) LAI, fl, 0 ct , r^, * w , z Q , 7 “ 6 

Sensible heat (Win 2 ) LAI, 9, tf cr , W 9 wi z Q , r~6 

Min. surface temperature (K) *, C s , kC s 


«*>!■ 


no 


Table 4.8a. Parameters selected at the PUPS forest site 

based on a threshold of |3<r| 


. Outliers 

Evaporation (mm/yr) 

None 

Runoff (mm/yr) 

None 

Sensible heat (Wm' 2 ) 

None 

Min. surface temperature (K) 



Table 4.8b. Parameters selected at the PILPS forest site 

based on a threshold of |2<r| 


Outliers 

Evaporation (mm/yr) r^, 6, LAI, $ a 

Runoff (mrn/yr) r^, o, LAI. * W) e CT 

Sensible heat (Wm* 2 ) r^, 6, LAI. * w , e a 

Min. surface temperature (K) *, C s , ^C s 


Ill 


Table 4.9. Identified important parameters (ranked fiom left to right) at 
the PILPS grassland site for the annual quantity metrics 



Primary 

Secondary 

Evaporation (mm/yr) 

LAI, 0, 0 C j., r^^, 8^ 

z 0’ W 

Runoff (mm/yr) 

LAI, 0, 0 CT , r^-^, 0 W 

z 0’ W 

Sensible heat (Wm'^) 

LAl, 0, k, 0 cr , r^-^, 0^ 

z 0’ r min d 

Min. surface temperature (K) 

K, Cg 

rtCg 


Table 4.10. Identified important parameters (ranked from left to right) at 
the PILPS forest site for the annual quantity metrics 



Primary 

Secondary 

Evaporation (mm/yr) 

r min’ * 

LAI, 0 W , 0 cr 

Runoff (mm/yr) 

r min’ S 

LAI, 0v7> Ar 

Sensible heat (Wm*^) 

r min’ ^ 

LAI, 8yf, 0 cr 

Min. surface temperature (K) 

K, Cg 

«Cg 
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Table 4.11. Results for b-, D s , and W s from one-at-a-time 
sensitivity analysis at the PUPS grassland site 


TfK) SE(Wm 2 ) ETfmm/yr) Rfmm/vr) 


Control run 

282.12 

47.87 

-5.10 

617.21 

648.83 

bj a 0.01 

282.11 

48.04 

-5.24 

619.44 

646.59 

bj = 0.03 

282.11 

48.00 

-5,20 

. 618.84 

647.19 

bj = 0.3 

282,12 

47.49 

-4.77 

612.20 

645.42 

b- a 0.3 

282.13 

47.20 

-4.53 

608.47 

652.84 

D 8 a 0.01 

282.12 . 

47.87 

-5.10 • . 

617.21 . 

648.84 

Dg = 0.03 

232.12 

47.87 

-5.10 

617.21 

648,73 

Dg = 0.1 

282.12 

47.87 

-5.10 

617.21 

648.83 

Dg = 0.3- 

282.12 . 

47.87 

-5.10 

617.21 

648.84 

tl 

O 

282.12 

47.87 

-5.10 . 

617.21 

648.83 

Wg a 0.3 

232.12 

47.S7 

-5.10 

617.21 

. 648.83 

Wg a 0.5 

282.12 

47.87 

-5.10 

617.21 

648,84 

Wg = 0.7 

282.12 

47.87 

-5.10 

617.21 

648.84 
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Table 4.12. Results for b-, D s , and W g from one-at-a-time 
sensitivity analysis at the PUPS forest site 


T(K) ^LwEfWm 2 ) S5fWm~ 2 ) ETfmm/v^ Rfmm/vr) 


Control run 

300.29 

104.55 

25.93 

1357,56 .. 

1909.50 

jr 

II 

o 

o 

300.23 

104.88 • 

25.62 

1361.88 

1905,20 

f-4 

6 

ii 

* •-» 

300.29 

103.81 

26.64 

1347.72 

1912.17 

bj = 0.3 

300.31 

102.17 

28*17 

1326,32 

1939.46 

bj = 0.5 

300.32 

101.00 .. 

29.28 

1310.96 

1955.85 

D s = 0.01 

300.29 

104.55 . 

25.93 

1357.57 

1909.52 

D s . — 0.05 , 

300.29 

104.55 

25.93 

1357.59 

1908.76 

D s = 0.1 

300.29 

104.55 . 

25.93 

1357.59 

1909.49 

D s = 0.3 

300.29 

104.55 

25.93 

1357.59 

1909.52 

W s = 0.1 

300.29 

104.55 

■ 25.93 . 

1357.59. 

1909.50 

CO 

o 

II 

g 

300.29 

104.55 

25.93 

1357.59 

. 1909,52 

W s = 0.5 

300.29 

104.55 

25.93 

1357.59 

1909.52 

$ 

II 

O 

300.29 

104.55 

25,93 

1357.59 

1909.51 
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Table 4.13. Results for and *cr from one-at-a-time sensitivity 

analysis at the PUPS grassland site 


T(K) ^UEfWnT 2 ) SHfWm' 2 ) ETfmm/yr) Rfmm/vr) 


r mia 

= 170 (s/m) 

282.08 

49.33 

-6.33 

635,94 

630.08 

r min 

= 230 (s/m) 

282.14 

46.61 

-4.0i 

600.93 

665,11 

0 CX : 

= 0.306 

282.07 

49.40 

-6.33 

636.94 

629.09 

*CT = 

= 0.383 

282.14 

46.77 

-4.15 

602.94 

663.06 


Table 4.14. Results for r m - n and 0 CZ from one-at-a-time sensitivity 

analysis at the.pJ!_r_- T..rc3t site 


TfK) Pw UEfWin 2 ) SHfWin' 2 ) ETfmm/yr) Rfmm/vr) 


r mia ~ 

: 140 (s/m) 

300.28 

105.72 

24.83 

i372.65 

1894.45 

^mia = 

' 160 (s/m) 

300.30 

103.78 

27.02 

1342.42 

1924.62 

o 

II 

0.33 

300.29 

104.69 

25.79 

1359.42 

1907.67 

cF 

h 

1! 

0.39 

300.29 

104.22 

26.24 

1353.36 

1913.73 
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Table 4.15a. Results of 32 runs for the sum of the monthly absolute 
difference metrics at PILPS grassland site 


Runs 

DET (mm/mo) 

DSH (Wm-2) 

DR (mm/ mo) 

i 

16.63 

64.24 

305.98 

2 

20,61 

43.76 

150.76 


37.08 

63.64 

85.12 

4 

25.40 

56.90 

241.58 

5 

72.84 

103.72 

304.62 

6 

51.24 

89.62 

255.21 . 

7 

221.08 

216.26 

370.96 . 

Q 

193.93 

198.58 . 

203.31 

9 

45.94 

63.29 

222,62 

10 

60.06 

97.01 

- 86.56 

11 

203.59 

134.95 

281.13 

12 

181.72 

108.56 

181.71 

13 

86.06 

42.06 

228.90 

14 

52.01 

27.95 

193.04 

15 • 

88.08 

32.85 

339.31 

16. 

75.86 

44.12. 

316.58 

17 

159.53 

133.31 

164.16 

18 

196.28 

151. ii 

244.78 

19 

33.58 

48.43 

272.18 

20 

94.90 

57.14 

231.74 - 

21 

• 105.99 

64.09 

224.48 

22 

145.01 

90.45 - 

—202.57 

22 

32.23 

79,87 

172,13 

2.4 

25.34 

46.68 

299.28 

25 

235.15 

212.29 

265.03 

26 - 

294.23 

256.96 

337.92 

27 

74.25 

77.93 

174.52 

28 

93.83 

109.25 

248.26 

29 

97.34 

99.12 

'200.27 

30 

117.00 

118.92 - 

180.77 

31 

16.75 

46.84 

46.39 

12 

20.69- 

64 . 19 

219.76 
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Table 4.15b. Results of 32 runs for the sum of the monthly absolute 
difference metrics at PILPS forest site 


Runs 

DET{ mm/mo) 

DSH(Wm~2) 

DR (rim/mo) 

1 

123.55 

98.36 

—799.03 

2 

196.60 

150.74 

361.73 

■5 

148.74 

102.65 

427.95 

4. 

115.90 

88.49 

545.81 

5 

208.12 

197.59 

544.46 

6 

129.81 

125.41 

1049.71 


404.26 . 

363.98 

723,21 

3 

359.80 

328.41 

836.55 

q 

254.42 

235.86 

496.20 

io 

244.72 

218.28 

294,07 

11 

309.24 

259.61 

488.83 

12 

134.73 . 

157.56 

477.93 

13* 

139.22 

111.73 

1023.77 

14 

■ 171.24 

139.37 

432.10 . 

15 

124.53 

95.31 

1269.52 

16 

133.40 

113.36 

559,87 

17 

194.30 

152.18 

920.04 

13 

252.99 

200.27 

432.92 

19 

189.14 

143.04 

1006.10 

20 

243.64. 

199.74 

465.76 

21 

86.25 

78.25 . 

568.41 

22. 

117.43 

92.22 

463,68 

23 

235.04 

220.51 

479.92 

24 

125.21 

125.69 

794.89 

25 

300.71 

271.90 

620.20 

26 

333.78 

340.14 

1276.18 

27 

143.39 

131.39 

365.44 

23 

134 . 13 

173.34. 

955.33 

29 

113.06 

106.89 

706.26 

30 

153.94 

143.47 

381.22 

31 

91.77 

95.93 

295.19 

32 

155.27 

135.91 

524.92 
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Table 4.16a. Parameter effects for the sum of the monthly absolute 
difference metrics at PUPS grassland site 


index Mo. 

DET 

(mm/mo) 

DR 

(mm/ mo) 

DSH 

(Wm-2 

1 

4.36 

- 4.00 

4.90 

2 

- 13.33 

7.26 

- 16.98 

. ■ ! 

■J - 

- 25.93 

16.47 

- 19.59 

4 

17.20 

- 12.-38 

1.73 

5 

21.52 

- 17.70 

16.82 

6 

4.10 

0.55 

1.79 

7 - 

- 20.31 

_ 3.24 

- 16.31 

3 

25.28 

- 19.95 

26.13 

9 

- 1.39 

21.09 

- 0.69 

10 

0.96 

14.58 

1.88 

11 

- 9 . 5.1 

- 58.60 

- 2.40 

12 

- 11.10 

29.06 

- 6.81 

13 

- 3.43 

1.93 

- 5.43 

1.4 

1.30 

4.30 

9.81 

15 

17.95 

59.74 

11.71 

' 16 

10.30 

14.96 

28.66 

17 

5.48 

20.52 

10.26 

13 

- 3.50 

- 28.32 

0.01 

19 

- 52.77 

- 25.56 

- 53.43 

110 

- 3.^0 

- 4.42 

37.52 

111 

1.46 

19.34 

6 .H 

24 

- 11.42 

4.30 

- 20.38 

25 

- 95.29 

- 26.73 

- 57.51 

2 . 3 ,* 

45.37 

11.25 

19.34 

35 

- 57.95 

- 65.59 

- 34.94 

56 

- 3 . 63 . 

57.22 

- 6.73 
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Table 4,16b. Parameter effects for the sum of the monthly absolute 
difference metrics at PILPS forest site 


Par Index No. DET 

(ram/mo) 

DR 

(mm/mo) 

DSH 

(Wm-2) 

1 

4.81 

-55,11 

4.17 

2 

4.88 

-9.55 ' 

4.49 

3 

-45.10 

45.01 

-28. 13 

4 

-2.07 

-15,82 

3.88 

5 

-16.77 

-4.65 

-11.02 

6 

10.59 

-10.15 

5.61 

7 

-23.32 

14.30 

-20.62 

8 

47.51 

-81.28 

45.64 

9 

31.90 

-287.82 

28.92 

10 

12.15 

236.33 

10,42 

11 - 

-5.6/ 

-123.50 

-6.35 . 

12 

-•23.43 

68.23 

-15.41 

13 

-12.45 

-15.86 

-12.46 

14 

- 10.81 

9.64 

9.93 

15 

26.79 

96.78 

22.17 

16 

58.27 

48.86 

56.04 

17 

10.23 

49,11 

7.46 

18 

-10.35 

-34.54 

-9,92 

19 

-70.24 

-17.67 

-77.64 

110 

13.58 

14.95 

19,51 

111 

11.96 

13.26 

8.87 

24 

-59/20 

-27.07 

-55.14 

25 

-34.24 

-50.62 

-24.46 

23 . 

42.81 

8.75 

33.54 

25 

-56.66 

-273.44 

-48.52 

56 

-18.33 

146.33 

-16.53 
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Table 4.17; Identified important parameters (ranked from left to right) 
at the PILPS grassland site for the sum of the 
monthly absolute difference metrics 


. ; Importance 

Evaporation (mm/mo) LAl0 cr .TO. ^nua LAI t LAI0 W * ifl&cx) 

Runoff (mm/mo) {Ml, r^LAI), W s , b.W s (*E^) 

Sensible heat (Wm' 2 ) LAI$ cr (32£), LAI$ W (33^), r m ^LAI (5^) 


Table 4.18. Identified important parameters (ranked from left, to right) 

at the PILPS forest site for the sum of the 

monthly absolute difference metrics 


Importance 

Evaporation (mm/mo) 
Runoff (mm/mo) 
Sensible beat (Wm'') 

bj, b^ (Ml, r^tAl), D s , b^TO, Wj 

LaIS, v (9®cr)> f min^ (®CT*)> f min^ 
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Table 4.19. Results of model effects at the PILPS-grassland site 



TfK) 

PwL^fWm* 2 ) 

SHfWm" 2 ) 

ETfmm/vr) Rfmm/vr) 

Std. Case 

282.12 

47.87 

-5.10 

617.21 

648.83 

Case (1) 

282.18 

43.80 

0^02 

565.51 

700.54 . 

Case (2) 

232.64 

48.30 

-7.90 

623.72 

642.32 . 

Case (3) 

280.35 

113.64 

-57.11 

1459.56 

— — — 

Case (4) 

280.51 

137.02 

-83.00 

1761.69 

— 

Table 4.20. 

Results of model effects at the PILPS forest site 


TfK) 

PwWSfWm' 2 ) 

SHfWm* 2 ) 

ETfmm/yr) 


Std. Case 

300.29 

104.55 

25.93 

1357.59 

1909.50 

Case (1) 

300.34 

90.60 

41.84 

1176.49 

2090.55 

Case (2) 

300.50 

157.04 

-25.31 

2034.11 

1232.97 

Case (3) 

299.78 „ 

146.09 

-9.65 

1900.05 

— 

Case (4) 

297.11 

484.35 

-332.33 

6293.91 

— 


Precipitation (mm/month) Precipitation (mm/month) 

0 200 ; 600 ' 0 200 400 
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Pig, 4.3 Parameter effects at the PILPS forest site. 
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Fig. 4.5 Effects of parameter bj on runoff at the PILPS grassland site. 
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Fig. 4.8 Parameter effects at the PILPS grassland site with the absolute 
difference as metrics. 
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Fig. 4.Q Parameter effects at the PILPS forest site with the absolute 
difference as metrics. . 





CHAPTER 5 A PARAMETERIZATION OF SPATIAL VARIABILITY OF 
PRECIPITATION IN THE TWO-LAYERYIC MODEL 


In Chapters 2-4, the structure o£ a two-layer VIC model was developed, 
and its performance and parameter sensitivities were evaluated. In the two 
applications, the model was driven with spatially-aVerage precipitation which 
was justifiable, since both sites were small. The FIFE site is a 15x15 km^' 
region and the ABRACOS site is essentially a point. For a large area (e.g., a 
GCM grid cell), however, the effects of subgrid scale spatial variations of 
precipitation on surface energy fluxes, Soil moisture, and runoff production may 
be significant. In practice, the subgrid land surface variations have been largely 
ignored in GCM land surface schemes. Most GCMs, for instance, assume 
uniform soil characteristics within a grid cell, and ignore spatial variability in 
precipitation. The two-layer VIC model described in Chapter 2 is a simple 
approach for representing subgrid variability in soil properties. In this chapter, 
an extension of the model to incorporate a representation of subgrid variability in 
precipitation is described. With this new representation, the effects of subgrid 
scale spatial variability of precipitation On surface fluxes, soil moisture, and 
runoff can be examined. The results of this method are compared with those 
obtained using an exhaustive pixel-based application of the two-layer VIC model, 
and those obtained by applying uniform spatial average precipitation to the two- 
layer VIC model. 

5.1. Introduction 

The feedbacks from land surfaces to atmospheric general circulation 
models (GCMs) are important determinants of regional and global climate. 
Land surface schemes used in the GCMs have to deal with spatial variations that 
occur at scales smaller than a GCM grid cell, in the same way that the 
atmospheric part of the GCM must parameterize subgrid scale atmospheric 
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variability. Natural land surfaces are heterogeneous within the scales resolvable 
by CCMs. There are two major types of heterogeneities. One is related to the 
subgrid scale hydrologic and topographic heterogeneity of the land surface itself; 
the other is related to the subgrid Scale Variability of meteorological inputs, Such 
as Precipitation, downward shortwave and long- wave radiation* wind speed, 
temperature and humidity. Among the meteorological inputs, the subgrid scale 
variability in precipitation is particularly important (Raupach 1993). Blyth et 
d. (1993) showed that correct prediction of the total evaporation can be obtained 
by using simple averages of surface parameters when comparing the results from 
a one-dimensional model With those from a three-dimensional mesoscale model. 
However, their results indicated that a reasonable partitioning of the total 
evaporation into, transpiration, evaporation from intercepted water, and 

evaporation of the bare soil cannot be obtained unless the spatial distribution of 
rainfall is considered at least. 

The spatial variability in precipitation has been widely recognized to have 
a major effect on evaporation, soil moisture variability* and runoff production 
(for example, Warrilow et al. 1986* Shuttleworth 1988, Entekhabi and 
Eagleson 1989, Famiglietti and Wood 1990, Pitman et al. 1990, Henderson- 
Sellers and Pitman .1992). Two approaches could be taken, to incorporate spatial 
variability in precipitation in a model. One is the pixel-based approach which 
discretizes precipitation over a Spatial domain. The. work of Famiglietti. et .al. 
(1992), and Wigmosta et al. (1994), among others, for example, falls into the 
- pixel-based category. Although such a pixel-based representation is able to 
account for the spatial variability in precipitation throughout a grid cell (or a 
catchment) in a straightforward manner, the associated computation time and 

data demands using this method make it untouchable_for implementations 
within ..GCMs. 

Another option is the statistical-dynamic approach to representing the 
spatial variability in precipitation. The advantage of this approach is that, if an 
appropriate statistical distribution is assumed, it can result in a closed form 
solution which would be computationally much less demanding than pixel-based 
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approaches. Also, this approach requires less data, since only the statistical 
distribution of precipitation needs to be specified, and not the values associated 
with’ specific spatial locations. Previous applications Of the statistical-dynamical 
approach axe briefly reviewed below. 

Warrilow et al. (1986) combined a subgrid precipitation distribution with 
a constant maximum surface infiltration rate -to estimate runoff from a grid celL 
They assumed that over a fraction, p, of a grid cell (or an area), the point 
precipitation intensity, Pj, is exponentially distributed and can be expressed as, 

frPi) = -p^exp( — prj-) > 0 < p < 1 (5.1) 

where P m is the grid cell average precipitation generated in the GCM. The 
runoff from the grid cell is then given by 


q. = 


*-J(P i -F*)ltPi)dP i '=P m «p(~ 

F* 



(5.2) 


where F* is the maximum surface infiltration rate. 

Shuttiewdrth (1988) derived an expression for canopy, throughfall based on 
the assumption that the precipitation rate over a fraction p of a grid cell was 
expressed by Eq. (5.1). Assuming that C m is the difference between the canopy 
storage capacity and the Water stored on the canopy divided by the model time 
step, the throughfall over the grid cell is then given.by 


P t = ,i.J (Pi-C m )f(Pi) dP; = 
Cm 


Pm exp( - 


pC 

~T 


m 


m 


(5.3) 


and the runoff over the grid cell becomes 
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q-P m 


eX p[— ^ 


(F*+C m ) 
^m J 


(5.4) 


Entekhabi and Eagleson (1989) represented the subgrid hydrologic 
processes m a GCM land surface scheme by assuming the same precipitation 
distribution (Eq. (5.1)) and combined it with a two-parameter gamma 
probability density function (pdf) of the surface layer point effective relative soil 
saturation to describe the Spatial heterogeneity in surface soil moisture. By 
assuming independence of the point precipitation intensity, Pj, and the surface 
layer point effective relative soil saturation, s, they derived a general relationship, 
for runoff rate (q) for the entire GCM grid during a time step as, 


/fiq = 


^ oo oo 

I (Pi -f )f(Pi)dPj f s (s)ds + f P i f(P i )dP i fg(s)ds (5.5) 

f* 10 




where f is the infiltrability (infiltration rate) of the first soil layer, f(Pj) is 
defined in Eq. (5.1), and f$(s) is a two-parameter, gamma pdf of s. 

Famiglietti and Wood (1991) considered the subgrid , scale variability in 
topography, soils, soil, moisture and precipitation by combining the same point 
precipitation distribution (Eq. (5.1)) with the distribution of the topography-soils 
indices from Topmodel (Beven and Kirkby 1979). They obtained an expression 
of the expected value of the depth of infiltration excess runoff for a large area 

E [Q inf ]> as. 


E I Q iofl=| } (Pi-f*)f(Pi)dPif 2 ( 2 )d 2 (5.6a) 

0 «?. 

where 2*=ln*(bTg/To tan ^). *i is the infiltration rate, 2=ln(bTg/TQtan/3) is 
the local value of the topography-soils index of Topmodel, 2 * represents the 
critical value of the Topmodel index at which saturation occurs, and f z (2) is the 
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pdf of 2 , which is assumed to follow a three-parameter gamma distribution 
according to Sivapalan et al. (1987) and Wolock et al. (1989). The expected 
value of the depth of saturation excess runoff* E[Q .], is expressed as, 

Saw 

oo oo z* oo 

E t Q satl = J J WiW + } } (Pi - Si)f(Pi)dPi f z ( z )dz (5.6b) 
z‘ 0 OS; 

where Sj is the storage deficit. The first integral in Eq. (5.6b) represents runoff 
generated on those areas that are saturated at the start of a time step, the 
second integral represents the runoff, generated on those areas that become 
saturated during a time step. 

Pitman et al. (1990) incorporated the parameterization suggested by 
Warriiow et al. (1986) and Shuttleworth (1988) into the BiosphereTAtmosphere 
Transfer Scheme (BATS) of Dickinson et al. (1986) to study the sensitivity of 
evaporation and runoff due to Eq. (5.1) within a grid cell. When incorporating 
Eq. (5.1) in BATS, they assumed that the soil moisture and intercepted water 
were distributed uniformly within a grid cell at. the end of a time. step. They 
found out, on the basis of results from a sensitivity study, that the monthly 
distributions of evaporation . and runoff were . quite sensitive to the spatial 
precipitation distribution and the fractional coverage ft of precipitation. 

Among the parameteriZations discussed above, only the ones by 
Entekhabi and Eagleson (1989), and by Famiglietti and Wood (1991) considered 
the effects of interactions of subgrid spatial variabilities in precipitation and 
subgrid variabilities in land surface characteristics. However, the runoff 
computed by both of their models (Eq. (5.5) or Eq. (5.6)) is a point average over 
the fraction ft of the grid cell on which precipitation falls. In other words, the 
runoff given by these models is equivalent to the runoff that would be generated 
by assuming that each point within the fraction ft is independent from each other 
and has the same statistical distribution of precipitation, soil moisture, and/or 
topography-soil index. Thomas and Henderson-Sellers (1991) pointed out that 
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these treatments ignore the absolute spatial location . within the grid cell of the 
precipitation. 

In this chapter, a derived distribution approach which combines the 
spatial subgrid variability in precipitation with subgrid variability in other land 
surface features, including soil moisture capacity, is described* To distinguish 
the statistical-dynamic approach used by Entekhabi and Eagleson (1989), and 
by Famiglietti and Wood (1991) from the one described here, we will refer their 
methods as point statistical-dynamic approach. The model developed here is 
described , as a one-dimensional statistical-dynamic model by using a derived 
distribution approach. 


5.2. One- di mens ional statistical-dynamic 

Assume that within a grid cell (or an area), precipitation, infiltration 
capacity, and/or other features (or attributes) only vary along one direction, 
which is arbitrarily taken as the x axis, and they axe kept constants along its 
orthogonal (y) axis (one-dimensional concept), where x and y axe scaled to give 
x-y — area. Ideally then, the general relationship for .runoff rate Q ^ of 
the area with a precipitation coverage n during a tim e step would be, 



[P(x) — f (x)] -f(Px)dPx f(z x )dz x + 


P x,4 z x,4 


| | P ( x ) *f(Px)dPx f(z x )dz x } *dx = p- 

P x,3 Z x,3 


dQ^, 0 < x < 1 (5.7) 


where P(x) and f (x) are precipitation and infiltration capacity within the 
fraction of area ^ in which rain occurs, f(P x ) and f(z x ), both of which vary with 
x, are pdfs of P x =P(x) and z x respectively, and z x is a variable which varies 
along x axis. This variable could be an effective relative soil saturation like that 
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in Eq. (5.5), or a topography-soils indices like that in Eq. (5.6), or some other 
variable, depending on the specific formulation. P x2 , z x 2 , P* 4 , and ± x 4 are 
the upper limits of the integrals in Eq. (5.7), while P x i, ’ z x v P x3 , and z x 3 
are the lower limits of the integrals in the same equation. It can be Seen that 
Eq. (5.7) reduces to a point statistical-dynamic approach (i.e.* Eqs. ( 5 . 5 ) and/or 
(5.6)) if the integrands in the bracket {} do not vary along the x axis. Since it is 
difficult to determine f(Px) and f(zx), and also to evaluate the triple integrals, 
we will simplify Eq. (5.7) Without loss of the one-dimensional treatment of 
subgnd spatial variabilities! In the derivations below, four assumptions are 
made: 


( 1 ) . The precipitation P is a one- dimensional function varying along the x axis 
(i.e., P(x)) within the. fractional coverage / 1 . This assumption is, in a sense, 
equivalent to assuming that storms are distributed as circles around the storm 
centers. . 

(2) . At the end of each time step t (i.e., at the beginning of the next time step 
t+ 1 ),. the soil, moisture content of each strip ”ydx” (defined in Section 5.2.1) 
within the fractional coverage ft of the Same, vegetation cover becomes the same. 
This assumption avoids the necessity for tracking the storm center movement. 

(3) . Prior to the beginning of the next storm, the soil moistures . over the 
fractional coverage n are assumed to be the same as the moistures over the non- 
rainfall fractional coverage of 1 — p, which is . accomplished by spatial averaging. 
This assumption becomes reasonable in practice if the inter-arrival time between 
two storms is long enough so that the recently wetted soil drains to a comparable 
moisture level to that which was not covered by the storm. Here we define the 
inter-arrival time as the time between two storms whose magnitudes are above a 
specified threshold (taken to be 1 mm/hr in this study). If a storm with 
magnitude below a specified threshold occurs, the soil moistures over /i and 1 - p 
are not averaged. It should be noted that if storms smaller than the specified 
threshold occur, they are not ignored, but from the standpoint of the soil 
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moisture distribution, they are treated as a Continuation of the previous storm. 
The effect of this assumption is that the track of previous storm centers relative 
to the present storm need not be Specified. 


(4). Each strip of area ydx is assumed to have an identical infiltration capacity 
distribution, which is defined by Eq. (2.20). 

5.2.1. One-di mftnsi onal representation for bare soil 

As discussed previously, Eq. (5.7) has to be simplified to obtain a one- 
dimensional runoff representation that could be implemented into GCMs. In the 
following, we use the direct runoff concept of the VIC model to obtain an 
expression for dQ^ d of Eq. (5.7). Let us begin with the simplest , case by 
assuming that only bare soil is. present in a grid cell (or an area) with a fraction . 
coverage of precipitation p at a time step t. If we discretize the area within the 
fractional coverage /i into infinitesimal strips of area ydx (see Fig. 5.1), then the 
precipitation rate, within each such strip is a constant. From .the fourth 
assumption and Eqs. (2.25a) and (2.25b), the . direct runoff dQ d [N+l] 
(following the notation in Chapter 2, the land cover class N+l is defined^ as bare 
soil) from Strip ydx due to precipitation P(x) can be expressed, as 

dC W N+1 l = ( Pl M _ “aT + “4t — “ ) dx » io+ p ( x ) • At > i m (5.8a) 


. 3- II- ,<->■“ n 


i 0 +P(x)-At < i m (5.8b) 

where, as defined in Chapter 2, Wj is the maximum Soil moisture content of 
layer 1, i Q and i m are the infiltration capacity and maximum infiltration 
capacity respectively, b- is the infiltration shape parameter, and At is the time 
step; and W^fN+l] is the soil moisture content in layer 1 within /i. Therefore, 


138 


the toted direct runoff. at time step t from p is, 


f W? W JN+l] 

V N+1 i = * ■ • J m*> ■ + 




(5.9) 


J m 


where the first term represents the direct runoff generated when ig + P(x) - At > 
im is satisfied, and the second term represents the runoff when ig + P(x) • At < 
im is satisfied. The integral limit, a, represents the location, x, where ig + 
P(x) • At = im if ig + P(l) • At < i m . If ig + P(l) • At > i m , then a = 1 in Eq. 
(5.9). 

Q 

In Eq. (5.9), Wj, bj, i m and At are constant. Based on assumption 2, 
the soil moisture W^[N+1] in Eq. (5.9) does not vary along the x axis. In 
addition, from Eq. (2.20) we can obtain the following relationship between ig 
and W pl [N+l], 


*0 


im[l~(l~ 


1 

W^lN+l]^ 1 +% i 
“Wf } 1 


(5.10) 


Thus, ig is independent of x. Therefore, only P(x) varies with x in Eq. (5.9). 
If the expression of P(x) is known, then the direct runoff Q^N+l] with the 
effects of spatial precipitation and infiltration can be calculated. For a constant 
precipitation tate within p, Eq. (5.9) reduces to Eq. (2.25). 

For evaporation, it is reasonable to ignore the effects of spatial variations 
in precipitation within the fractional coverage p for big storms, which is 
equivalent to assuming that evaporation is small during storm periods compared 
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with dry periods. For small storms, the effects of . spatial variations in 
precipitation are Small, and so can be ignored also. Therefore, the evaporation 
from bare soil (E^) can be determined by Eq. (2.22) multiplied by p, as if the 
precipitation were uniformly distributed over p. Tnat is 


E pi = /i : E i- 


(5.11) 


The drainage term Q^IN+I] is determined by Eq. (2.2?) multiplied by p: 

Q y l 2 [N+H » y°- . (5.12) 


The water balance in layer 1 over the fraction ft where precipitation occurs is 
then, 


w 1 [N+l]=W; 1 [N+l]+(p 


P(x)dx^Q^j[N+l]-Q^j 2 [N+l]-E p i)-At ( 5 . 13 ) 


where W^N+l] and W^j[N+l] are the soil moisture content-at the end and the 
beginning of each time step in layer 1 within p respectively. S imi larly the 
subsurface runoff Q^fN+l] and the water balance for the lower layer within p 
can be expressed as, 


Qpb( N+1 l = ^ S S W$ W p2l N+1 )> 0 ^ W^[N+1] * w s W 2 (5.14a) 

W , 2 [ N + 1 | + ( D m ^)( ^ 1 ^ W 2 )2)t 


r s ”’2 


WJj--W s W§ 


W^ 2 [N+ 1 ]>W s W§ (5.14b) 
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and 


w J[N+l]=W^[N+l]+(Q pl2 [N+ll * Q^iN+1] - E„ 2 ) ■ it (5.15) 

where W^N+l] and W^fN+l] are the soil moisture Content at the end and the 
beginning of each time step in layer 2 within n, respectively. Other terms have 
the same definitions as used in Chapter 2. 

For the area within the dry (no rain) fraction 1 - n, the evaporation from 
the bare soil drainage from layer 1 to layer 2, Q^.^^N+l], 

subsurface runoff Q^ ^^[N-M], and soil moistures in layer 1 and layer 2, 
W(i- M )i[N + i], ^[N-fi], are calculated in the same way as described 

above, except using W^.^jjN+l], W (i. /i )2l N+1 ]’ 311(1 ( 1- ^) instead of 

W^i(N+l], W^2[N+1]j and /i in.Eqs. (5.11) -(5.15). The direct runoff in this 
case, is zer0 * 

According to assumption 3, the soil moisture within area p and 1 - p is 
averaged immediately prior to the beginning of next storm that is larger than the 
specified threshold so that the entire area has the same, soil moisture content 
when next, storm arrives.. Thus, the soil moisture in layer j (j=l, 2) 
immediately prior to the above threshold follow-up storm is, 


W;,j[N+l] - ^ • (W^[N+l)+W^. rij [N+l)), (5.16a) 

W (Mi [N+11 = (‘-/‘)-(W^lS+l]+wf 1 .„) j [N+l)). (5.16b) 


5.2.2. One- dimensional representation for vegetation cover 

For the area with different vegetation covers, the approach described 
above applies except that the rainfall rate P(x) is replaced by the throughfall 
rate. It is assumed that the throughfall rate P t (x) is equal to P(x) minus 
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interception Where interception is calculated by Eq. (2.3). The approach for 
incorporating different vegetation covers and bare soil within n is the same as 
described in Chapter 2. 

Like the bare soil case, the soil moistures for different cover classes over 
the wet and dry fractions n and l-/u are averaged (or "smeared”) prior to the 
beginning of the next storm that is larger than the specified threshold. That is, 
the soil moistures of bare soil (W^N+1], j=l, 2) and each vegetation type n 

( W /ij[ n ]’ n=1 > 2 > * • *> N) over the fractional coverage of precipitation /a, and 
the soil moistures of each such cover type over the dry fractional coverage 1 
are averaged, so that the entire area (or grid cell) has the same Soil moisture 
content at the onset of the next storm. For the fractional area a covered by 
precipitation, the average soil moisture V/,j (j=l, 2) can be expressed as, 

N+l , 

w pj - E C v [n].W p j[n]. (5.1.7a) 

n=l 

Similarly, the average soil moisture over the dry fractional area 1 

can be expressed as, 

N+l . 

w (l-/i)j- E c v[a]*W(i_ /J ) j [n]. (5.17b) 

n=l 


The average soU moisture of layer j (j=l, 2) for the entire area, at the beginning 
of the next storm is, 

W/ijM * M '( W pj +W (l- P )j)> n = 1, 2, • • ., N+l (5.18a) 

W (l^)j fnJ = (1 " ^ ‘ (W /ij +W (l-/i)j)’ “ = 1, 2, , N+l. (5.18b) 

There are two places in which the soil moisture is averaged. The first is over 
strips of area ydx within n at each time step for each cover class using Eq. (5.13). 
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This avoids tracking the movement of a storm center during a storm. This 
averaging results in equal soil moisture within the same surface cover type within 
ft. It is different among the different surface cover classes, and it is different 
from the soil moisture for the dry fraction 1 — ft for the Same surface cover class. 
The second averaging, described by Eq. (5.16) or Eqs. (5.17) and (5.18), is 
earned out only at the beginning of the next above-threshold storm. 

The two advantages of this derived distribution approach axe that it 
avoids the need to identify the specific area that receives -a given precipitation 
rate within the fraction ft of an area (or a grid cell), and it considers the spatial 
variability of precipitation within the fractional area covered by the storm. This 
is- because for P(x) (or P t (x)) within a. strip ydx, P(x) (or P t (x)) is considered to 
be a constant, and the spatial variability of the infiltration function (Eq. (2.20)) 
is considered over a strip area ydx with the same initial soil moisture. The dry 
fraction l- ft of the area (or grid, cell) is taken to. have no precipitation 
throughout, the storm. The fraction ft is assumed to be a constant within each 
storm, but it can vary from stofm to storm. 

5.2.3. Derivation of precipitation function along x axis 

The one-dimensional statistical-dynamic . model requires that the 
precipitation function P(x) be known. In this section, an appropriate form of 
P(x) is derived. 

Eagleson (1984), Eagleson and Wang (1985), and Eagleson et al. (1987) 
have reported studies on the fraction ft of a grid area that is affected by 
precipitation reaching the surface. _ Warrilow et al. (1986) assumed that over a 
fraction of a gnd cell (of an. area), the point precipitation intensity is 
exponentially distributed (Eq. (5.1)). This exponential distribution assumption 
for precipitation seems appropriate in some cases as shown by an analysis of 
hourly observed precipitation data by Abdulla (1987). Abdulla (personal 
communication 1993) also found that precipitation over a large area in Oklahoma 
appears to follow the exponential distribution. Using the continental-scale radar 
data over northwestern Europe, Collier (1993) found that for frontal' 
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precipitation, the exponential relationship (Eq. (5.1)) may not be appropriate. 
Instead, the log-normal relationship used by Foufoula-GeOrgiou and Wilson 
(1990) might be more appropriate. Gao and Sorooshian (1994) investigated ten 
years of hourly precipitation data over three GCM grid squares covering the 
southeast, southwest, and north central U.S. each having a size of 8 degrees of . 
latitude by 10 degrees of longitude. They found that the exponential distribution 
assumption may be questionable and that the statistical .patterns of precipitation 
depend on the locations and the interactions between atmospheric conditions and 
various land surface characteristics. On the basis of the previous work, two 
possible forms of the precipitation function P(x) (based on exponential and log- 
normal distribution) will be considered. 

First, following Warrilow et al. (1986), the precipitation intensity is 
assumed to be described by Eq. (5.1) within p. The percentage (x) that receives 
a precipitation rate larger than or equal to precipitation rate P over the .fraction 
/i of a grid cell (or an area) can then be expressed as (see Fig. 5.1a), 


.. ^ area that receives precipitation rate > P 
* “total area of„ ^(P) 


P 

= 1- }f(P,)dP; 
0 





(5.19) 


The inverse of Eq. (5. 19.) is then, 

p (x) = -- ^P-ln(x) , 0 < x < 1. (5.20) 

In Eq. (5.20), P(l) = 0, that is, the probability of the area fraction p that can 
receive precipitation greater than or equal to zero is one. This result is due to 
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the-form of Eq. (5.1). 

For the second case where the log-nOrmal distribution is assumed (which 
is arguably more appropriate for frontal precipitation), the precipitation 
intensity function is expressed as (Collier 1993), 

f ( p i) ~ [ ^(1 + Pj)]^ } (5.21) 

where all the symbols have the same meaning as before* By approximating the 
integral of the probability distribution* Collier (1993) obtained, 


F(P) = [f(Pi)dPj = 1- exp{--J!-[ Ml + P)] 2 }• 

6 m 


(5.22) 


Thus similarly to Eq. (5.19), the percentage (x) that receives, a precipitation 
rate larger than or equal _to. precipitation rate P over the fraction n can then be 
expressed as, 


_ area that receives precipitation rate > P _ . 

total area of p — l-r(P) 

= exp{ -- Ji- [ ln(l + P)] 2 } . (5>23) 

* m 


The inverse of Eq. (5.23) is then, 


1 

P(x) = exp { - ln( x )]"^"} - 1, 


0 < x < 1. 


(5.24) 


Here again Eq. (5.24) yields P(l) = 0. 

For the exponential distribution case, by substituting Eq. (5.20) into Eq. 
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(5.9), we obtain, 


/*• 




r P W? 

= /i*| { -L + 


a.isg!!,.*. 


/*• 


, Pmw., W{ W^N+1] W 

1 I* m W- aF+ At +AT1 1 i~ J }dx 


— Pmr f* (-*r + 


Wf w 4 [n+i] 


At ' At 


) + ^j, 1 .!oz few^ , 1+ V 


(5.25) 


The integral in Eq. (5.25), when 0 < a < 1, can be computed as, 


[1- 


iQ — ^-ln(x)-At i+^i 


— Om) 


_ p 
P m 


■(l+bj) ‘“^o 


] dx = 


r »m - *0 + — jp" ln(x) • At , 1 + b i , 

l J <ix 


a 


J m 


d l+ h i F^At Ct-ila+io) 


1 PnvAt* * e 


tl 


dt 


0 


m ' At .dt 


(5.26) 
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where t « i m - i^+ ln( x ). A t, M d tl = i m -i 0 + h(a). ^ {acti 

tl = i m -io + — ln(a) is equal to . zero, since a is determined by 

i 0 + P ( x ) = *in ^ ^Qf"F(l) * At < i m (see Section 5.2.1). 

An analytical expression can be obtained for tbe case where bj is an 
integer in Eq. (5.26). Otherwise, numerical integration can be used to calculate 
the integral in Eq. (5.26), or to expand the term expfe-ii— t) into a power 
senes which converges Within It I < oo: 


exp(i 


n-t 


m- 


At 


) = 1 + 


fi-t 


+ *[• 


1 r M-t -,2 


m * At 2! Tm • At 


r+* 


i 1 r in , 

+ 5PF*^tl +' 

Itl < oo . 


(5.27) 


For the log-normal distribution case, by Substituting Eq. (5.24) into Eq. 
(5.9), Q^IN+l] can be obtained through numerical integration. 

5.2.4. Estimation of. fractional coverage of precipitation 

As for the fraction p, WarriloW et al. (1986) noted that p may be on the 
order of 0.95 and 0.60 for large-scale and convective rainfall, respectively. 
Currently the U.K. Meteorological Office (UKMO) Sets p to be 1.0 for large-scale 
rainfall, and. 0.3 for convective rainfall in their GCM. The observed spatial 
variability of total storm depth for air. mass thunderstorm rainfall in Arizona 
(Eagleson et al. 1987) supports a Wetted fraction of 0.5 to 0.66. While WarriloW 
et al. (1986), Entekhabi. and Eagleson (1989), Pitman et al. (1990), and 
Famiglietti and Wood (1991) all used a prescribed fraction p of a grid cell area, 
Eltahir and Bras (1993) found significant temporal variability in the fractional 
coverage of rainfall. They presented .a procedure for estimating the fraction p as, 

** = E(FV) (5-28) 

Where E(Pp) is the areal mean precipitation over the rain-covered fraction p_of 
the grid cell. They suggested that P m can be estimated by 
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Pm = ~a X' a\ ■ aT (5.29) 

where V is the volume of rainfaU simulated by the chmate model within a grid 
cell area, AX, and AT are the spatial and temporal resolutions of the model 
respectively. By invoking the ergodicity assumption, they estimated E(P M ) by 
the mean of the rainfall rate at a point using rainfall records from a single 
location. Although the volume of rainfall, V, could be taken from a numerical 
weather prediction model, this can result, in large biases; it is probably more 
realistic to use an average of observed station data instead of Eq. (5.29) to 
estimate P m . Based on the observed hourly rainfall data in the Southeast .of 
France, Braud et al. (1993) found that the fractional area covered by rainfall 
exceeding a fixed threshold is highly correlated with the mean areal rainfall rate. 
Thus, this suggested that n may also be estimated from the radar rainfall data. 

5.3. Testing of the derived distribution approach 

As described previously, there are four major assumptions involved in 
deriving the one-dimensional statistical-dynamical model. In this section, two 
computer experiments are designed to test these, assumptions by applying Eq. 

(5.20) (i.e., the exponential . distribution) to simulate spatially varying 
precipitation. The first is a "brute force” experiment (Figs. 5.2a and 5.2b) where 
a grid cell (or an area) is divided into LxM pixels. In the tests that were 
performed, there were LxM=250Q pixels (subgrid elements). At each specific 
time, the average precipitation depth of an area (or grid cell) defined the mean 
of the exponential precipitation, distribution. Thus, the exponential precipitation 
distribution over the area at each time was determined. The precipitation rates 
for the 2500 pixels Were taken randomly from the exponential distribution (Eq. 

(5.20) ) and were then assigned to each pixel. In generating the rainfall field, the 
variable x in Eq. (5.20) was obtained from a uniform random generator. Wit hin 
each pixel, the precipitation rate was assumed to be the same. The procedure of 
assigning the precipitation rate to each pixel Was repeated at each time when 
there was rainfall. Using this random method, it is possible for a pixel to have 
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large value of precipitation in one time period and be followed by a random 
quantity drawn from the exponential distribution in the next time period. No 
attempt was made to incorporate the spatial correlation in precipitation 
amounts. It was assumed that one infiltration capacity distribution characterized 
the entire area. The infiltration capacity distribution for each pixel Was obtained 
by randomly sampling from the infiltration capacity distribution for the entire 
area (Fig. 5.2b). These pixel hydrologic properties were kept unchanged during 
the simulations for a fixed spatial precipitation coverage y. Since the 
precipitation rate was assumed to be the same within each pixel, the two-layer 
VIC model described in Chapter 2 Can be applied. By applying the model 
described in Chapter 2 on a pixel' by pixel basis, the direct runoff, evaporation, 
sensible heat flux, surface temperature and soil .moisture were calculated for 
each pixel and then were aggregated to obtain results for the entire area. 

In the second experiment (Fig. 5.2c), the derived distribution approach 
described in. Section 5.2. was conducted to compute the direct runoff, 
evaporation, sensible heat flux, surface temperature, and soil moisture for the 
same, area as in Experiment 1. In this experiment, the precipitation over the 
area y followed an exponential distribution with the same distribution 
parameters used in Experiment 1. 

Finally, the results of the direct runoff, soil moisture content, 
evaporation, sensible heat flux, and surface temperature from Experiments 1 
and 2 were compared with the results obtained by assuming P(x) = constant 
over the area. 

The same forcing data and model parameters as were used in the Chapter 
4 sensitivity analysis are used here, except for the vegetation cover parameter . 
(C v ). which is taken as 1.0 for both the grassland and forest sites. Other 
parameter values are as listed in Table 4.1. Since the model does not represent 
variations in spatial snow properties, only summer month simulations for the 
grassland were conducted. At the forest site, snow does not occur, so the 
simulations were conducted for an entire year. For simplicity, the precipitation 
fractional coverage parameter y was constant in each simulation, but three 
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different values of n (0.3, 0.6, and 1.0) were investigated. In the following 
discussion of these experiments, we use ’‘pixel-based'.’, "derived distribution”, 
and "average” to represent Experiment 1 (the "brute force" experiment), 
experiment 2, and the one with P(x) = constant within the area, respectively. 
All simulations were conducted at an hourly time Step. 

At the forest site, the 12^month hourly results for the latent and sensible 
heat fluxes, surface temperature, soil moisture of the upper zone, and soil 
moisture for both zones of the three experiments are shown in Figs. 5.3-5.14 for. 
n = 0.3. All experiments started with the same initial soil moisture which was 
taken as 50% of the maximum soil water content. On the same figure, the 
hourly area average precipitation time series is also plotted for comparison. Six- 
month (Feb., Mar., April, Oct., Nov., and Dec.) hourly runoff time series for 
months with maximum runoff peaks higher than 0.6 mm/hr are shown in Fig. 
5.15 for /i=0.3. 

In January (Fig. 5.3), many small storms with intensities less than 3.0 
mm/hr were pres. nt. The evaporation simulated by the. derived distribution 
approach is very close to that simulated by the pixel-based approach, while the 
average method predicts much larger evaporation. The reason is that the 
average approach assumes that precipitation is intercepted within the entire area 
while only 30% of the area lias intercepted precipitation. in the pixel-based and. 
derived distribution approaches. Therefore, the average approach results in 
much more interception evaporation than the other two approaches, especially 
when the precipitation rate is small. Likewise, for the sensible heat flux and 
surface temperature, the derived distribution approach gives results that are 
much closer to the pixel-based approach than the average approach. The ratio of 
the monthly sum of the absolute difference between the average and pixel-based 
approaches to the monthly sum of the absolute difference between the derived 
distribution and pixel-based approaches is shown in Table 5.1a for the latent and 
sensible heat fluxes, and surface temperature. From this s umma ry table, it can 
be seen that all the ratios are greater than 3, indicating that the difference 
between the derived distribution approach and the pixel-based approach is much 
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smaller than that between the average and pixel-based approaches. 

From Fig. 5.3, it is also seen that the derived distribution approach 
better approximates the soil moisture content of the tipper zone and also the 
moisture of the upper and lower zones combined than the average approach. 
The monthly maximum relative errors of the upper zone moisture content are 
6.8% and 7.7% for the derived distribution and average approaches respectively, 
and 0.4% and 1.3% for the combined upper and lower zone moisture content for 
the derived distribution and average approaches respectively (Table 5.3a). 

In January, the largest runoff peak is small (0.04 mm/hr) due to the ' 
small precipitation rate and the. low soil moisture content. In Table 5.2a, a. 
summary ratio for runoff is. calculated in the same Way as for the latent and 
sensible heat fluxes, and surface temperature shown, in Table. 5.1. For months . 
with the maximum daily runoff rate less than 1.0 mm/day, the. ratio is 
uninformative and is not calculated. In Table 5.2a, the monthly 
hourly peak runoff is shown for the pixel-based, derived distribution, and 
average approaches as well. 

In February (Fig. 5.4), there were not as many storms as in January, 
but their magnitudes were larger. Even though the evaporation from the average 
approach is still larger than the pixel-based and derived distribution approaches, 
the difference is not as large as in January due to smaller interception. The 
sensible heat flux and surface temperature also varied less among the three 
approaches than in January. However, the total runoff from the average 
approach Was significantly undersimulated compared with the pixel-based one 
(Fig. 5.15). During the two big storms at the end of the month, the average soil 
moisture over the grid area is higher in the average approach, but the overall 
saturated area is smaller than in the pixel-based and derived distribution 
approaches, and thus the average approach generates much less direct runoff. 
The derived distribution approach produced simulations of runoff (Fig. 5.15), 
evaporation, sensible heat, surface temperature, and soil moisture (Fig. 5.4) 
that matched the pixel-based, results quite closely. The ly.tios V>r maoff, latent 
and sensible heat fluxes, and surface temperature are 7.35 (Table 5.t’ii), 1.59, . 
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1.88, and 1.56 respectively (Table 5.1a). The monthly maximum relative errors 
of the upper zone moisture content axe 7.3% and 19.1% for the derived 
distribution and average approaches respectively, and 0.8% and 4.5% for the 
combined upper and lower zone moisture content for the derived distribution and 
average approaches respectively (Table 5.3a). 

The simulated results for March (Fig. 5.5) are quite similar to the results 
for February. At the beginning of the month, both the pixel-based and derived 
distribution approaches predict much larger runoff (Fig. 5.15) than in January 
even though the storm magnitudes were similar to those for January. This is 
because of the wet soil from the large storms at the end of February. The 
average approach again gives much smaller runoff in this case. In comparison, 
the storm with peak rate more than 3 mm/hr occurring on the 19th of the month 
results in very small runoff in all the three approaches (only about one fifth of . 
the runoff, that occurred at the beginning of the month) due to the loss of soil, 
moisture by evaporation during the period. The ratios for runoff, latent and 
sensible heat fluxes, and surface temperature are 2.41 (Table 5.2a), 1.72, 1.69, 
and 1.47. respectively (Table 5.1a). The monthly maximum relative errors of the 
upper zone moisture content are 8.0% and 14.0% for the derived distribution and 
average approaches respectively, and 0.7% and 4.9% for the combined upper and 
lower zone moisture content for the derived distribution and average approaches 
•respectively (Table 5.3a). 

In April, storms .were smaller than in February, but larger than in 
January and March. Again, the derived distribution approach approximates the 
runoff (Fig. 5.15), surface fluxes, and soil moisture (Fig. 5.6) much better than 
the average approach. The ratios for the latent and sensible heat fluxes, and 
surface temperature are all greater than 2.0 (Table 5.1a), and greater than 4.0 
for runoff (Table 5.2a). Although the evaporation simulated by the average 
approach is always larger than those, from the pixel-based and derived 
distribution approaches, the much smaller runoff from the average approach 
results in higher soil moisture than for the other two approaches. The monthly 
maximum relative error of the upper zone moisture -content is greater than 20% 
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for the average approach, but less than 10% for the derived distribution 
approach (Table 5.3a). For the combined moisture content, the monthly 
maximum relative errors are 0.7% and 6.2% for the derived distribution and 
average approaches respectively (Table 5.3a). 

May, June, July, and August (Figs. 5.7 — 5.10) is the dry period. There 
is little runoff during these months, and the largest monthly runoff peak is less 
than 0.5 mm/day, except in June where it is 4.32 mm/day. During this period, 
the average approach continues to simulate larger evaporation due to more 
moisture, except for August. The soil moisture content for the average 
approach decreases and. approaches the soil moisture from the pixel-based 
approach. During this time, the derived distribution approach still gives, better, 
approximations than the average approach with all the ratios for the latent and 
sensible heat fluxes, and surface temperature greater than 1.0 (Table 5.1a), 
except for August. For soil moisture, the comparison (Table 5.3a) shows that 
the derived distribution approach results in smaller monthly maximum relative 
errors than the average approach for both the upper layer and combined layers, 
except for the upper layer soil moisture in June and July, and the combined soil 
moisture in August, where the average approach gives slightly smaller relative 
errors. 

In August (Fig. 5.10), there is almost no precipitation and the soil is 
quite dry. The soil moisture from the derived distribution . approaches closer to 
the pixel-based approach than the average approach. The evaporation from the 
average approach is overestimated during the first half of the month, and also 
during the three very small precipitation events due to larger interception. The 
evaporation from the derived distribution approach is also overestimated in this 
month. . 

The reason that the derived distribution approach overpredicts the latent 
heat flux might be attributable to the way of calculating the soil moisture stress 
factor gsmW (Eq. 2.14). In the pixel-based approach, the subgrid spatial 
variability of soil moisture is included in the calculation of the soil moisture 
stress factor gsmMi but it is not in the derived distribution and average 
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approaches. When the soil is dry, the pixel-based approach may result in a 
number, of large gsm[ I1 ] ,s over a large portion of the area due to, for example, 
less soil moisture than the wilting point or critical point, while the, derived 
distribution and average approaches give only one smaller lumped g S m[ n ] value. 
Therefore, both the derived distribution and average approaches predict larger 
evaporation than the pixel-based approach. This subgrid spatial effect would be 
more significant when the soil is drier since for example, the g sm [n] would.be 
the same (i.e., 1.0) for the pixel-based, derived distribution, and average 

approaches if the entire area is saturated. 

Although both the derived distribution and average approaches use the 
same lumped expression for g sm [n], the average approach simulates smaller 
evaporation than the derived distribution approach during this month (Fig. 
5.10). This is because the upper zone soil moisture is drier in the. average 
approach than in the derived distribution approach. Thus, the average 

approach can result in a larger soil moisture stress factor gsmM than the derived 
distribution .approach and so simulates smaller evaporation. Therefore, the 
average approach simulates evaporation better than the derived distribution 
approach but for the wrong reasons. 

From September to December (Figs. 5.11-5.14), the soil becomes wet 
again as a result of more storms, some of which are large. Therefore, the effect 
of subgrid spatial variability on g sm [n] is reduced among the three approaches. 
The evaporation from the average approach is oversimulated during most of 
these months and the runoff (Fig. 5.15) is significantly undersimulated for the - 
large precipitation events compared with the pixel-based results. The derived 
distribution approach, however, gives much better approximations. All the 
ratios for the latent and senaible heat fluxes, and surface temperature are 
greater than 1.0 (Table 5.1a) during these months, and the ratio for runoff is as 
high as 5.68 in October (Table 5.2a). The monthly maximum relative error for 
the upper zone soil moisture is about. 10% or less for the derived distribution 
approach, but about 20% in some months for the average approach (Table 5.3a). 
Similarly, the derived distribution approach results in monthly raaidmum 
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relative errors for combined soil moisture smaller than the average approach, 
except for September. 

For the. case ju^O.6, ' similar results are obtained as for /i=0.3, except that 
the underestimation of runoff and overestimation of evaporation by the average 
approach decreases, for the larger /i value. The ratios for the latent and sensible 
heat fluxes, and surface temperature are greater than 1.0 for 9 out of 12 months, 
and are all above 0.9 for March (Table 5.1b). The ratios for runoff are greater 
them 1.0 except for March and November (Table 5.2b). The monthly maximum 
relative errors for the upper zone soil moisture are about 10% or less for the 
derived distribution approach, while they are as high as 21.3% for the average 
approach (Table 5.3b). The monthly maximum relative errors for the combined 
soil moisture are small fcr both approaches; the largest relative errors over the 
year axe 2.5% and 3.2% for the derived distribution and average approaches 
respectively. 

At the grassland site, which has a drier, climate them the forest site, 
simulations for (i— 0.3, 0.6, and 1.0 were conducted for the summer months 
(May -September), with the initial soil .moisture equal W the maximum soil 
moisture content for ^=0.6, and half of the maximum soil moisture content in 
the other two cases. 

For ••■=0.3, Figs. 5.16 and 5.17 show that May and June receive less 
precipitation than the remaining three months, and the maximum hourly runoff 
is 0.01 mm/hr and 0.03 mm/hr for May and June, respectively (Table 5.5a). 
The average approach .ovei simulates evaporation during the precipitation event 
due to large interception. The derived distribution approach produces much 
better approximations for evaporation, sensible heat flux and surface 
temperature. The ratios f nr the latent and sensible heat fluxes, and surface 
temperature are all greater than. 3.0, and are as high as 10.53 (Table 5.4a). The 
derived distribution approach also results in smaller monthly maximum relative 
errors for the upper zone and combined soil moisture (Table 5.6a). 

Evaporation is significantly overestimated by the average approach in 
July and August (Figs. 5.18 and 5.19) during the storm period, and thus the 
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sensible heat flux and surface temperature are underestimated; while the derived 
distribution approach yields the evaporation, sensible heat, and surface 
temperature similar to the pixel-based approach, with ratios greater than 5.0 
(Table 5.4a). Although the simulated runoff is small during this period with a 
maximum rate of 0.09 mm/hr and 0.22 mm/hr for July and August respectively 
(Table 5.5a), both the derived distribution and average approaches give larger 
runoff thaii the estimates from the pixel-based approach. This is because both 
the derived distribution and average approaches always have some saturated area 
when the soil moisture in the upper zone is greater than Zero according to the 
variable infiltration capacity formulation for the upper layer. Therefore, when 
there is some small rainfall, the precipitation that falls onto the saturated area 
will generate runoff in both cases. For the pixel-based approach, however, it is 
not necessary that such a saturated area always exists, and thus there is no 
runoff from light precipitation unless the antecedent soil moisture is saturated. 
When the precipitation rate is moderate or large, the effect of the perpetually 
saturated area in the derived distribution and average approaches becomes 
negligible. The derived distribution approach again results in smaller monthly 
maximum relative errors for the upper zone and combined soil moisture (Table 
5.6a). Figure 5.20 shows the results for September in which the derived 
distribution approach again produces better simulations than the average 
approach (Tables. 5.4a-5.6a). 

For the case p=1.0, similar comparison results are obtained for the three 
approaches as for the /*= 0.3 case, except that underestimation of runoff and 
overestimation of evaporation by the average approach decreases as ft increases, 
as was the case for the forest site. The ratios for the latent and sensible heat 
fluxes, surfr.ce temperature, and runoff of the five m iths (May - September) 
are all greater- than 1.0, except for June where the ratio of latent heat flux is 
0.97 (Tables 5.4b and 5.5b). Although for most months, the monthly maximum 
relative errors for the upper zone and combined soil moisture are slightly greater 
for the derived distribution approach than for the average approach, all the 
values of the relative errors are small, with the largest less than 3% (Table 
5.6b). 
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To investigate the effect of initial soil moisture on the surface fluxes and 
runoff, for the drier climate, a comparison of the three approaches with initial 
soil moisture set to the maximum for p=0.6 Was made. Comparisons (Tables 
5.4 -5.6) show that the results are similar except that the runoff from the pixel- 
based approach is smoother than from the other two approaches in the p=0.6 
case (plots not included). This is because more water goes into lower zone and 
comes out as subsurface runoff in the pixel-based approach. Tables 5.4c and 5.5c 
show. that all the ratios for the latent and sensible heat fluxes, surface 
temperature, and runoff are greater than 1.0, except for the ratio of runoff in 
June which is 0.95. Also, the average approach shows, larger relative errors for 
the soil moisture in general (Table 5.6c). 

Tables 5.7 and 5.8 show the annual average latent and sensible, heat fluxes 
at the forest site and the five-month average latent and sensible heat fluxes at 
the grassland site, respectively, for. the three approaches. For all of the cases 
with /i=0.3 and 0.6 at the forest site, and ^=0.3, 0.6, and 1.0 at the grassland 
site, the derived distribution approach gives better approximations to the pixel- 
based approach than the average approach. At the forest site, the average 
approach gives about 20 Wm' 2 bias for the latent and sensible heat fluxes, and 
the derived distribution approach gives about 5 Wm 2 bias. At the grassland 
site, the average approach gives about 20 Wm* 2 bias for the latent heat flux and . 
over 35. Wm*- 2 bias for the sensible heat flux fo- ^=0.3. The derived distribution 
approach, however, gives less than 5 Wm* 2 difference for both latent and 
sensible heat fluxes for ^=0.3. 

Although the experiments conducted here use hypothetical forcing data 
with specified initial soil moisture, the differences for the surface fluxes, surface 
temperature, runoff, and soil moisture between the average and pixel-based 
approaches indicate that the effect of spatial subgrid scale variability in 
precipitation can be significant. Comparison of the three approaches shows. that 
the derived ._ distribution approach approximates the pixel-based method 
reasonably well in terms of surface fluxes, surface temperature, runoff, and soil 
moisture. 


157 


5.4.._S tih nu ary of the denYed distribution approach 

This chapter has described a one- dimensional statistical-dynamic model 
with a derived distribution approach that accounts for the effect of subgrid 
spatial variability of rainfall. The approach Was tested against a pixel-based 
approach, and an approach which applied spatial average precipitation 
uniformly over the area. The comparison of the three approaches for the PILPS 
forest and grassland sites shows that: 

(1) . The derived distribution approach approximates the pixel-based approach 
much better than the average approach in the, simulations of surface fluxes, 
surface temperature, runoff^_and soil moisture; 

(2) . The effect of the spatial subgrid scale variability in precipitation is 
significant. The inclusion of the spatial variation of precipitation results in less 
evaporation, especially for small storms,, due to less interception of rainfall, 
and more runoff, especially for moderate and large storms, due to a larger 
portion of saturated area. The soil moisture is lower in general when the spatial 
variation in precipitation is considered for moderate or large storms because 
more runoff occurs, while the soi 1 moisture is higher for small storms due to less 
evaporation. 


(3) . The ratios of the monthly sum of the absolute difference between the 

average and pixel-based approaches to the monthly sum of the . absolute 
difference between the derived distribution and pixel-based approaches for the 
latent and sensible heat fluxes, surface temperature, and runoff are greater than 
1.0 for most cases. The monthly maximum relative errors for the upper zone soil 
moisture . and for the combined soil moisture are smaller for the derived 
distribution approach than for the average approach for most cues. 

(4) . Comparisons with the pixel-based approach show that the annual average 
and five-month average latent and sensible heat fluxes at the forest and grassland 
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sites for the average approach have biases of about 20 Wm* 2 , while the derived 
distribution approach has biases of about 5 Wm* 2 . . The average approach 

resulted in over 35 Win. 2 bias in the sensible heat flux for fi= 0.3 at the grassland 
site. 

(5). The four assumptions described in Section 5.2 incorporated in the derived 
distribution approach seem justifiable based on the test results. Although the 
third assumption, which averages the soil moisture over the fractional coverage 
(i and the non-rainfall fractional coverage 1-p, may seem -questionable, the 
results Suggest that it. is defensible. Under a number of common situations, this 
assumption may not be critical. One is for the case of many small storms. In 
this case most of. the rainfall is intercepted by the vegetation canopy. Thus, the 
soil moisture over the fractional coverage n would, remain close to that over the 
non-rainfall fractional coverage 1^. The second case is when several large 
storms occur, with most of the area covered by at least one storm. In this case, 
the soil moisture will be close to saturation everywhere- The third case is where, 
the soils are relatively well drained, which will tend to minimize spatial 
differences in soil moisture. 


Since the derived distribution approach is conceptually simple and has a 
closed form, it is computationally tractable. It appears to be a practical way of 
representing subgrid scale precipitation variation in a GCM land surface 
parameterization. 


Table 5.1. Comparison of surface fluxes and temperature between the derived 
distribution ami average approaches for the PILPS forest site 



(T 6 ), and R e _ -g-, R„ _ -g-, Rq^_ -g- . 

^initial an ^ ^ represent the initial and maximum soil' moisture- of combined' upper and lower 
layers respectively. 



160 




<0 CT . 

5 - . >» 


5 2 
J. i 


>> >* 



WW 
ii ii 

V 

o 

2 


where hourly quantity is the total runoff Q, and Hq = -j~ • 

Qm(P). an< ^ Qim(^) represent monthly maximum hourly runoff (mm/hr) from the 

pixel-based, derived distribution* and> average approaches* respectively* 







Table 5.3. Comparison of soil moisture at the PILPS forest site 



W.’ — monthly maximum relative error obtained by — 1 : — ' x 100% 

YV^ixel 

where W* vg represents hourly soil moisture of layer 1 from the average approach. 

W' rr and W"* represent the monthly maximum relative error of soil moisture of combined 
upper and lower layers for the derived distribution and, average approaches respectively. 





Table 5.4. Comparison of surface fluxes anti temperature between the derived 
1 distribution and average approaches for the P1LPS grassland site 
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where hourly quantity can be evaporation (E), sensible heat flux (H), or surface temperature 
(T 6 ), and R e = R„, = -jp 11^= . 

W initial ant * wC rc P rcscnl the initial and maximum soil moisture of combined upper anti lower 
layers respectively. 


Table 5.5. Comparison of runoff at the PILPS grassland site 
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b = 53| hourly quantity from derived dist. appro, - hourly quantity from pixel-based appro. | 
where hourly quantity is the total runoff Q, and Rq = -jj- . 

Q m (P), Q in (D), and Q rn (A) represent monthly maximum hourly runoff (mm/hr) from the 
pixel-based, derived distribution, and average approaches respectively. 
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Table 5 JL — Comparison of annual average latent and sensible 
heat fluxes at the PILPS forest site 


/ i = 0.3 




latent 

sensible 

latent 

sensible 


(Wm* 2 ) 

(Win' 2 ) 

(Wm' 2 ) 

(Wm' 2 ) 

pixel-based 

66.5 

64.5 

73.7 

• 57.8 

derived dist. 

72.0 

59.3 

76.9 

54.7 

average 

87.7 

44.5 

87.7 

44.5 


Table 5.8. Comparison of five-month average latent and Sensible 
heat fluxes at the PILPS grassland site 



A «= 

= 0.3 


= 0.6 


= 1.0 


latent 

sensible 

latent 

sensible 

latent 

sensible 


(Wm* 2 ) 

(Wm* 2 ) 

(Wm 2 ) 

(Wm' 2 ) 

(Wm' 2 ) 

(Wm* 2 ) 

pixel-based 

50.8 

40.2 

73.3 

20.5 



derived dist. 

53.4 

39.7 

77.7 

18.5 


BUB 

average 

71.5 

2.0 

83.9 

- 8.6 

71.5 

2.0 


71.5 2.0 
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Fig. 5.1 Schematic representation of the derived distribution approach to 
represent spatial variability of rainfall for two-layer VIC model: (a) 
exponential precipitation distribution in x direction, and plan view 
of a grid cell (or area) with strips of ydx^ and ydx 2 * (b) two-layer 
VIC representation for the strips ydxj and ydxQ. 
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Fig, 5.4 Comparison of effective surface fluxes, temperature^ and soil moisture among pixel-based (solid), derived 
distribution (dashed); and average (dotted) approaches for February with /i=0.3 for PILPS forest site. 
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Fig. 5.5 Comparison of effective surface fluxes, temperature, and soil moisture among pixel-based (solid), derived 
distribution- (dashed) and average (dotted) approaches for March with /i=0.3 for PILPS forest site. 
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Fig. 5.6 Comparison of effective surface fluxes, temperature,, and soil moisture, among pixel-based (solid), derived 
distribution (dashed), and' average (dotted) approaches for April with /i=0.3 for PILPS forest site. 
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Fig. 5.9 Comparison of effective surface fluxes, temperature, and soil moisture among pixel-based (solid), derived 
distribution (dashed) and average (dotted) approaches for July with /i=0.3 for PILPS forest site. 






































Fig. 5.12 Comparison of effective surface fluxes,, temperature,, and soil moisture among pixel-based (solid), derived 
distribution (dashed) and average (dotted) approaches for October with /i=0.3 for PILPS forest site. 
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Fig. 5.14 Comparison of effective surface fluxes, temperature, and soil moisture among pixel based 1 (solid), derived 
distribution (dashed) and. average (dotted) approaches for December with ft- 0.3 for PILPS forest site. 
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Fig. 5.15 Comparison o£ runoff among pixel-based (solid), derived distribution (dashed) and average 
(dotted) approaches for /i=0.3 for the PILPS forest site. 
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Fig. 5.17 Comparison of effective surface fluxes, temperature, and' soil moisture among pixel-based (solid), derived 
distribution (dashed) and average (dotted) approaches for June with /i=0.3 for PILPS grassland site. 






2 51 - 



Fig. 5.18 G 
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Fig ? 5. ID Comparison of effective surface fluxes-, temperature, and soil- moisture among pixel-based (solid), derived 
distribution (dashed) and average (dotted} approaches for August with /<= 0.3 for PILPS grassland site. 
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CHAPTER 6 CONCLUSIONS AND RECOMMENDATIONS 

FOR FUTURE WORK 


A simple two-layer variable infiltration capacity (VIC-2L) model has been 
developed to model land-atmosphere . interactions at the scale of general 
circulation models used for numerical weather prediction and climate change 
studies. The model includes a canopy layer, which partitions the area of 
interest into N land surface cover types. For each land Cover type, the fraction 
of roots in the upper and lower layer is specified. Evaporation occurs via canopy 
evaporation, evaporation from bare soil (land cover class N+l), and 
transpiration, which is represented via a canopy and architectural resistance 
formulation. In the soil, the effects of spatial subgrid variability of soil moisture 
with hydrologically plausible runoff mechanism are represented through the 
upper layer and lower layer. The upper layer, which is designed to represent 
the dynamic behavior of the soil as it responds to r ainfall , is characterized by 
spatial distributions of infiltration and soil moisture capacities. The lower layer, 
which is used to characterize the seasonal soil moisture behavior, is spatially 
lumped and uses the Amo drainage representation. Drainage from the upper 
layer to the lower layer is assumed to be driven by gravity. The effect of the 
subgrid Spatial distribution of rainfall is accounted using a derived distribution 
approach. This approach results in a one-dimensional statistical-dyn ami c closure 
form based on four assumptions. The four assumptions are: (1) precipitation 
vanes along x axis. within a fractional coverage /*; (2) the soil moisture content 
of each Stnp within the fractional coverage fi of the same vegetation cover is 
averaged at the end of each time step; (3) the soil moisture is averaged over the 
fractional coverage n and the non-rainfall fractional coverage 1 - p prior to the 
beginning of next storm; and (4) each strip has an identical infiltration capacity 
distribution function. The VIC-2L model includes both atmospheric and 
hydrologic model parameters. 
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6.1. Conclusions 

. The VIC-2L model performed well for two applications where the climate 
regimes are quite different. The first application was to the FIFE site in central 
Kansas, where the climate was relatively dry, and the land cover was tall grass. 
The second site was in north central Brazil, which has a moist, tropical 
rainforest climate, and where the land cover is a ranch clearing surrounded by 
forest. 

At both the FIFE (35 days of data) and ABRACOS (59 days of data) 
sites, the VIC-2L model reproduced the latent heat flux well, with 80% and 
90%, respectively, of the daytime hourly peaks having relative errors less than 
15%. Although more than 80% of the days (FIFE and ABRACOS) had sensible 
heat flux daytime hourly peaks less than 200Wm*^, there were still about 60% 
and 55%, respectively, of the daytime hourly peaks with relative errors less 
than 15%. For the ground heat flux, the majority of days had daytime hourly 
peaks on the order of 50 Wm'^, and there were about 40% and 55% of the 
daytime peaks, respectively, with relative errors less than 15%. For the upper 
layer soil moisture, the comparison at the ABRACOS site (no observed soil 
moisture data were available at the FIFE site) showed that the largest relative 
difference between the model simulated and the observed, soil moisture was less 
than 8%. There were 70% of the days with surface temperature difference less 
than 2 °C at the FIFE site. For the streamflow at the FIFE site (not available 
at the ABRACOS site), the reproduction of dry period flows and the timing of 
the major peaks of the observed streamflow was satisfactory, although the 
magnitudes of the large peaks were subject to substantial errors. 

Since most rain over land is convective, the inclusion of spatial subgrid 
scale variability in precipitation should provide a better representation of land 
surface dynamics than the average precipitation representation for those areas. 
The effects of subgrid spatial variability of precipitation on surface fluxes, soil 
moisture, and runoff were incorporated into, the VIC-2L model. The model 
performance was evaluated by comparing hourly simulations of latent and 
sensible heat fluxes,, surface temperature, runoff, and soil moistures with the 
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pixel-based approach simulations for two different (hypothetical) climate 
regimes. In addition, the simulated latent and sensible heat fluxes, surface 
temperature, runoff, and soil moisture were compared with the results obtained 
by assuming constant precipitation over the area. The comparison showed that’ 
the derived distribution approach approximated the pixel-based approach quite 
well, and it Was superior to constant precipitation approach in terms of 
predicted surface fluxes, surface temperature, runoff, and soil moistures almost 
all the time. The simulations with constant precipitation over the area 
overestimated latent heat flux and underestimated sensible heat flux compared 
With the pixel-Wed simulations due to much mote interception evaporation, 
especially When the precipitation rate is small. Also, the constant precipitation 
experiment simulated much Smaller runoff, particularly for median and large 
precipitation, as compared, with the pixel- based results. In addition, the bias of 
the annual average and five-month average latent and sensible heat fluxes at the 
forest and grassland sites Was about 20 Wm* 2 from the average precipitation 
approach, while it was only about 5 Wm' 2 from the derived distribution 
approach. 

Finally, the VIC-2L model parameters were explored for two different 
climate regimes through both fractional factorial and one-at-a-time sensitivity 
analyses. Three metrics were chosen: annual total evaporation, annual total 
runoff, and annually averaged sensible heat flux. For these metrics, the leaf 
area index, porosity, and the minimum stomata! resistance were found to be 
important, while the soil thermal conductivity and soil heat capacity were found 
to be most sensitive for the minimum hourly surface temperature metric. These 
results are similar to those obtained by Henderson-Sellers (1992, 1993) using the 
BATS model. When the sum of the absolute difference of monthly evaporation 
and sensible heat flux between the control experiment and the 32 experiments- 
were used as metrics, the model Was found to be most sensitive to a similar set 
of parameters as those indicated by, using the annual total . evaporation and 
annually averaged sensible heat flux as metrics. However, for the absolute 
difference metrics, it Was the interactions of these parameters rather than the 
single parameter effects that were found to be most important. 
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The infiltration shape parameter, the fraction of. maximum subsurface 
flow, and the fraction of the lower layer maximum soil moisture were found to 
have minimal influence on the annual total runoff metric. Although the two 
fraction parameters redistribute the monthly runoff within a year, they did not 
change the total amount of the annual runoff. However, these three parameters 
and t heir interactions among themselves and with the parameters of porosity and 
wilting point were found to be sensitive when the sum of the absolute difference 
of monthly runoff from the control experiment was used as the metric. 


6.2. Recommendations for future work 

There are a number of research topics related to the current work that are 
suggested by this Work, and are Worthy of future exploration. Three suggestions 
for immediate follow up work related to representation of Subgrid scale 
variability of the soil moisture stress factor associated with transpiration, 
subgrid scale variability in snow properties, and model validation and 
implementation into GCMs are described briefly. 

6.2.1. Subgrid scale spatial variability in g r ,. 

From Chapter 2, it is seen that the representation of evaporation from 
bare soil (Eq. 2.22) accounts for subgrid spatial variability in soil moisture, but 
the effect of subgrid spatial Variability of soil moisture on transpiration from 
vegetation (Eqs. 2.12-2.14) is only partially accounted for. The weakness of Eq. 
2.14 is shown in Chapter 5 when comparing with the exhaustive experiment 
results for the forest site in August. Therefore, it is necessary to account for the 
subgrid spatial variability of Soil moisture on transpiration in a manner similar to 
that used for. bare soil evaporation. This can be done through the soil moisture 
stress factor g sm [n] (Eqs. 2.13 and 2.14) which depends on the water available in 
the root zone. Due to the definition of the upper layer and lower layer, the soil 
moisture in the upper layer is more spatially variable than that in the lower 
layer. Thus, the soil moisture stress factor gj^Jn] of layer 1 should have larger 
variations within a grid cell (or an area) than that of layer 2. We can, therefore, 


only consider the influence of subgrid Variability of soil moisture in layer 1 on 
ggujn]- A lumped soil moisture Content expression of g sm [n] can be used for 
layer 2 as before. From the definition of g sm [n] (Eqs. 2.14a-2.14c), we obtain, 
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dW 1 [n]-dw7 

dW^-dW^ 


1 

dA+| O dA 
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w cr 

where Wj_[n], , and Wj have the same meanings as in Chapter 2, and 

AcrM and- Ay^n] represent the fraction of the soil whose soil moisture is greater 
than or equal to the critical value above which the transpiration is not affected 
by the moisture stress in the soil, and the value associated With the permanent 
wilting point respectively, and (see Fig. 2.2), 


dWjJn] = ig-dA 
dW^ = 0 W • i • dA 


dWj* = 0 cr • i • dA 

where 0 W is the permanent wilting point (dimensionless), 0 cr is the critical value 
(dimensionless), ig represents the point infiltration capacity corresponding to the 
soil moisture Wj[n], and A is the fraction of the area for which the infiltration 
capacity is less than i (see Eq. 2.20). Let 

i w = 0 w -i , 
icr = ®cr'b 

then combining with Eq. 2.20, we have 


Thus, in Eq. (6.1), the second term can be expressed as* 
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The integrand of Eq. (6.4) can be either numerically integrated or expanded into 
a power series. The terms A^fn] and A cr [n] axe determined as follows. From 
Eq. 2.20, we can obtain the relationship between ig and W]Jn], 


iQ ~ im[l ~ (1 - 



(6.5) 


Equating Eq. (6.2) with Eq. (6.5), and Eq. (6.3) with Eq. (6.5) respectively, 
then, 

1 _ 

i0 = *w-im [l - (1 - A w [n]) 1 ] = i m [l - (1 — UijjiL) + 1 ] (6.6) 

Wj 

1 

i0 = *cfim[l-(l-A cr [n]) 7 ' ] = i m [l-(l-^L) + * ]. (6.7) 

WJ 


Therefore, if i m -0 w > ig> we have for A w [n], _ 


192 


1 -( 1 - 


Wj[n] 1+bj 


Aw[ n ] — 1 — [1 - ■ 


Wi 


7 w 


) 


] 


( 6 . 8 ) 


and if i m • 0 cr > ig, we have for A C f[n], 


A C r[n] — 1 - [1 -■ 


W-fnl 1+b i 

Wi , b i 

^cr . 


If im-^w < iQ’ thenAwW=l; while if i m • fl cr < i Q , A cr [n]=l. 


(6.9) 


6.2.2. Subgrid scale spatial variability in -a snowpack 

In . the midcontinents of North American and Eurasia, winter snow and 
spring jmelt are important components of the hydrologic cycle. Snowmelt in 
spring generally produces Spring peaks in streamflow. Numerous snowmelt 
models are described in the literature. These models can be grouped into three 
classes (Morris. 1985): regression models, lumped conceptual, models, and. 

distributed models. Each class of these models has its strengths and weaknesses. 
The subgrid spatial variability in snowpack properties has not been included into 
the present form of the VIC-2L model. However, in a context consistent with 
the approach described in. Chapter 5 for the subgrid spatial variability of 
precipitation, the subgrid spatial variability related to snowpack properties 
could be accounted for with. a. similar degree of complexity. 

One way to include the subgrid spatial variability of snow is to 
incorporate Donald’s (1992) findings of a log-normal Spatial distribution for snow 
depth and linear relationship of snow depletion curve into the snowmelt module 
(Wigmosta et al. 1994) of the VIC-2L model by using the derived distribution 
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approach described in Chapter 5. Based on snow Surveys conducted in southern 
Ontario, Canada, Donald found that snow depth can be described by a log- 
normal distribution within a shallow and often discontinuous SnOw covered area. 
In addition, he showed when Snowmelt is assumed to occur uniformly over the 
snowcovered area, the redistribution and accumulation of snow within an area 
classified by vegetation cover types can be summarized by linear snow depletion 
curves (SDC) for different cover types. The snow depletion curve, which can be 
characterized by two classes, forest and non-forest, describes the relationship 
between average snowcover depth and snowcovered area for a given vegetation 
cover type (Donald 1992). 

Two assumptions are needed to incorporate Donald’s findings. The first is 
that the area (or a grid cell) is relatively flat so that topographic effects such as 
slopes, aspects, and shading are insignificant over the snowcovered area (i.e., 
the net- radiation received over the snowcovered area. is relatively uniform). The 
second is that the rain falling on the snowcovered area will be redistributed in 
such a way that the heat advected to the snowpack by the rainfall will be . the 
same within the Snowcovered area. With these two assumptions, it can be 
assumed that the snowmelt will occur uniformly over the snowcovered area, and 
Donald’s linear relationship of SDC can then be applied. 

6.2.3. Validation and implementation into GCMs 

As seen in Chapters 3 and 5, the VIC-2L model was only applied to two 
specific sites, both of which were small. Therefore, the model should be applied 
to more sites with different soil type properties and climate conditions. In 
addition, observed data when available over a large area should be used, to test 
the one-dimensional statistical-dynamic model described in Chapter^, since the 
approach of including the subgrid spatial variability of precipitation (see Chapter 
5) was not tested using observed data. One possibility is to use the data from 
HAPEX-MOBILHY (Hydrological Atmospheric Pilot Experiment) which were 
collected over an area on the order of 100x100 kin. The observations Were 
conducted in 1986 with one Intensive Field Campaign (IFC) of duration several 
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months in southwestern France (Shuttleworth 1991b). The study area with 
multisite measurements consists of forest, agricultural land, and some wood 
cover area. The observed data include surface pressure, air temperature, 
dewpoint temperature, wind, vegetation information, radiation information, 
and streamflow (Goutorbe and Tarrieu 1991). 

To test the global performance of the VIC-2L model, the effects of 
feedbacks from GCMs should be investigated. Therefore, implementation of the 
VIC-2L model into a GCM is a logical next step in evaluation of the model. 
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